text:
# 1 Introduction

## 1 Introduction

Hadronisation, the transition of quarks into hadrons, is a strong interaction phenomenon which cannot yet be calculated from first principles within QCD. Monte Carlo event generators are used instead which rely on phenomenological models of this process. To some extent these models can be distinguished from each other by the shape of the predicted hadron energy distribution. Hadronisation of heavy quarks leads to a significantly harder hadron energy spectrum than for lighter quarks . Experimentally, heavy quark hadronisation is of special interest, because in this case the hadron containing the primary quark can easily be identified.

A precise measurement of the $`\mathrm{B}`$ hadron111All hadrons containing a $`\mathrm{b}`$ quark will be referred to as $`\mathrm{B}`$ hadrons throughout this paper. energy distribution allows the various hadronisation models available to be tested, and also helps to reduce one of the most important systematic uncertainties in many heavy flavour analyses. Earlier measurements of the $`\mathrm{B}`$ hadron energy distribution usually fell into one of three categories. 1) Some analyses were based on a measurement of the energy distribution of certain exclusive $`\mathrm{B}`$ hadron decays, mostly $`\mathrm{B}\mathrm{D}^{}\mathrm{}\nu `$, to constrain the $`\mathrm{B}`$ hadron energy as precisely as possible . However, this leads to small candidate samples and thus to a large statistical uncertainty. 2) Other analyses attempted to increase the sample size by extracting the energy distribution of leptons from inclusive $`\mathrm{B}\mathrm{}`$ decays. Unfortunately the modelling of the lepton energy spectrum introduces large additional systematic uncertainties . 3) The most precise results so far have been achieved through the fully inclusive reconstruction of the $`\mathrm{B}`$ hadron energy . The analysis presented here identifies $`\mathrm{B}`$ hadrons inclusively using secondary vertices.

## 2 Data sample and event selection

This analysis uses data taken at or near the $`\mathrm{Z}`$ resonance with the OPAL detector at LEP between 1992–2000. A detailed description of the OPAL detector can be found elsewhere . The most important components of the detector for this analysis are the silicon microvertex detector, the tracking chambers, and the electromagnetic calorimeter. The microvertex detector consisted of two layers of silicon strip detectors which provided high spatial resolution near the interaction region. The central jet chamber was optimised for good spatial resolution in the plane perpendicular to the beam axis222The OPAL coordinate system is defined as a right-handed Cartesian coordinate system, with the $`x`$-axis pointing in the plane of the LEP collider towards the centre of the ring and the $`z`$-axis along the electron beam direction.. The resolution along the beam direction was improved by the $`z`$ information delivered by the silicon microvertex detector (except in the first version present in 1992), by a vertex drift chamber between the silicon detector and the jet chamber, and by dedicated $`z`$-chambers surrounding the other tracking chambers. The central detector provided good double track resolution and precise determination of the momenta of charged particles by measuring the curvature of their trajectories in a magnetic field of $`0.435`$T. The solenoid was mounted outside the tracking chambers but inside the electromagnetic calorimeter, which consisted of approximately $`\mathrm{12\hspace{0.17em}000}`$ lead glass blocks. The electromagnetic calorimeter was surrounded by a hadronic calorimeter and muon detectors.

Hadronic events are selected as described in Ref. , giving a hadronic $`\mathrm{Z}`$ selection efficiency of $`(98.1\pm 0.5)\%`$ and a background of less than $`0.1\%`$. Only data that were taken with the silicon microvertex detector in operation are used for this analysis. A data sample of about 3.91 million hadronic events is selected. This includes 0.41 million events taken for detector calibration purposes during the years 1996–2000, when LEP was operating at higher energies.

A total of 23.81 million Monte Carlo simulated events are used: 16.81 million events were generated with the JETSET 7.4 generator , 2 million events were generated with HERWIG 5.9 , and 5 million events were produced by HERWIG 6.2 . The JETSET event sample includes 4.93 million $`\mathrm{b}\overline{\mathrm{b}}`$ events and 3.19 million $`\mathrm{c}\overline{\mathrm{c}}`$ events in dedicated heavy flavour samples. All other samples are mixed five flavour $`\mathrm{Z}\mathrm{q}\overline{\mathrm{q}}`$ event samples. The choice of important parameters of the event generators is described in . All Monte Carlo simulated events are passed through a detailed detector simulation . The same reconstruction algorithms as for data are applied to simulated events.

The analysis is performed separately for the data of different years, where detector upgrades, in particular of the silicon microvertex detector , and recalibrations lead to different conditions. Separate samples of JETSET Monte Carlo are available for all years. HERWIG Monte Carlo is only available for the largest homogeneous dataset taken in 1994, and therefore HERWIG-related studies are performed exclusively for this dataset.

In the 1993 and 1995 runs, part of the data was taken at centre-of-mass energies about $`1.8\mathrm{GeV}`$ above and below the peak of the $`\mathrm{Z}`$ resonance. The $`\mathrm{B}`$ hadron energy distribution is sensitive to energy losses due to initial state radiation prior to the annihilation process. Initial state radiation is heavily suppressed at and just below the $`\mathrm{Z}`$ resonance, but it has significant impact in the dataset taken at an energy of $`m_\mathrm{Z}+1.8\mathrm{GeV}`$. The latter samples are therefore treated separately, with Monte Carlo samples simulated for the appropriate energy, giving a total of eleven separate data and JETSET Monte Carlo samples.

## 3 Preselection of $`𝐙\mathbf{}𝐛\overline{𝐛}`$ events

The thrust axis is calculated for each event using tracks and electromagnetic clusters not associated with any tracks. To select events within the fiducial acceptance of the silicon microvertex detector and the barrel electromagnetic calorimeter, the thrust axis direction is required to satisfy $`\mathbf{|}\mathrm{𝐜𝐨𝐬}𝜽_𝑻\mathbf{|}\mathbf{<}\mathbf{0.8}`$, where $`𝜽_𝑻`$ is the thrust angle with respect to the beam direction.

To achieve optimal $`𝐛`$-tagging performance, each event is forced into a 2-jet topology using the Durham jet-finding scheme . In calculating the visible energies and momenta of the event and of individual jets, corrections are applied to prevent double counting of energy in the case of tracks with associated clusters . A $`𝐛`$-tagging algorithm is applied to each jet using three independent methods: lifetime tag, high $`𝒑_𝒕`$ lepton tag and jet-shape tag. This algorithm was developed for and used in the OPAL Higgs boson searches. A detailed description of the algorithm can be found in . Its applicability to events recorded at the $`𝐙`$ resonance peak was already shown in . The $`𝐛`$-tagging discriminants calculated for each of the jets in the event are combined to yield an event $`𝐛`$ likelihood $`𝓑_{\mathrm{𝐞𝐯𝐞𝐧𝐭}}`$. Both the jet b-tagging discriminant and $`𝓑_{\mathrm{𝐞𝐯𝐞𝐧𝐭}}`$ have values between zero and one and correspond approximately to the probability of a true $`𝐛`$ jet or $`𝐛\overline{𝐛}`$ event, respectively. For each event, $`𝓑_{\mathrm{𝐞𝐯𝐞𝐧𝐭}}\mathbf{>}\mathbf{0.2}`$ is required. The $`𝐙\mathbf{}𝐛\overline{𝐛}`$ event purity is 83% after this requirement, and the efficiency is 54% at this stage.

The $`𝐛`$ hemisphere tag efficiency obtained from Monte Carlo simulation is compared to the actual value in data using a double tag approach as described in . The efficiencies obtained this way in both simulation and real data are found to agree to within 5% in all subsamples. Nevertheless a correction is applied to the Monte Carlo efficiency to further improve the agreement.

## 4 Reconstruction of $`𝐁`$ hadron energy

The primary event vertex is reconstructed using the tracks in the event constrained to the average position of the $`𝐞^\mathbf{+}𝐞^{\mathbf{}}`$ collision point. For the $`𝐁`$ hadron reconstruction, tracks and electromagnetic calorimeter clusters with no associated track are combined into jets using a cone algorithm333Studies have shown that the cone jet-finder provides the best $`\mathrm{B}`$ hadron energy and direction resolution compared to other jet finders .  with a cone half-angle of $`\mathbf{0.65}\mathrm{𝐫𝐚𝐝}`$ and a minimum jet energy of $`\mathbf{5.0}\mathrm{𝐆𝐞𝐕}`$. The two most energetic jets of each event are assumed to contain the $`𝐁`$ hadrons. Only jets where the opposite hemisphere yields a $`𝐛`$-tagging discriminant of at least 0.8, corresponding to a $`𝐛`$ probability of about 80%, are used in the analysis. The distribution of the $`𝐛`$-tagging discriminant is shown in Figure 1.

Each remaining jet is searched for secondary vertices using a vertex reconstruction algorithm similar to that described in , making use of the tracking information in both the $`𝒓\mathit{\varphi }`$ and $`𝒓𝒛`$ planes where available. If a secondary vertex is found, the primary vertex is re-fitted excluding the tracks assigned to the secondary vertex. Secondary vertex candidates are accepted and called ‘good’ secondary vertices if they contain at least three tracks. If there is more than one good secondary vertex attached to a jet, the vertex with the largest number of significant444A track is called significant if its impact parameter significance with respect to the primary vertex is larger than 2.5. The impact parameter significance is defined as the impact parameter of a track divided by the uncertainty on this quantity. tracks is taken. If there are two or more such vertices, the secondary vertex with the larger separation significance with respect to the primary vertex is taken. Jets without an associated secondary vertex are rejected. This increases the $`𝐛`$ jet purity and improves the energy resolution of the $`𝐁`$ hadron reconstruction described in the following.

Weakly decaying $`𝐁`$ hadrons are reconstructed inclusively with a method described in an earlier publication . In each hemisphere defined by the positive axis of the jet found by the cone algorithm, a weight is assigned to each track and each cluster, where the weight corresponds to the probability that this track or cluster is a product of the $`𝐁`$ hadron decay. This weight is obtained from artificial neural networks  exploiting information from track impact parameters with respect to the primary and secondary vertices, and from kinematic quantities like the transverse momentum associated with a track or cluster, measured with respect to the cone jet axis. A list of all variables is shown in Table 1. The $`𝐁`$ hadron momentum is then reconstructed by summing the weighted momenta of the tracks and clusters. A beam energy constraint assuming a two-body decay of the $`𝐙`$ and the world average $`𝐁`$ meson mass of $`\mathbf{5.279}\mathrm{𝐆𝐞𝐕}\mathbf{/}𝒄^\mathrm{𝟐}`$  for the $`𝐁`$ hadron is applied to improve the energy resolution. The constraints lead to a biased energy reconstruction, particularly when the true $`𝐁`$ hadron energy is very small, as can be seen in Figure 2a. However, only a small fraction of the data sample is in the low-energy region affected by a large bias. For most events, in the peak of the $`𝐁`$ hadron energy distribution, the bias is small, and all biases are taken into account by the fitting procedures used in both the model-dependent and model-independent analyses. Possible systematic uncertainties arising from the biased energy reconstruction are discussed in Section 7 of this paper. The energy of the weakly decaying $`𝐁`$ hadron is expressed in terms of the scaled energy $`𝒙_𝐄\mathbf{=}𝑬_𝐁\mathbf{/}𝑬_{𝒃𝒆𝒂𝒎}`$, where $`𝑬_{𝒃𝒆𝒂𝒎}\mathbf{=}\sqrt{𝒔}\mathbf{/}\mathrm{𝟐}`$ is the LEP beam energy for the event. The quantity $`𝒙_𝐄`$ is restricted to values above $`\mathbf{5.279}\mathrm{𝐆𝐞𝐕}\mathbf{/}𝑬_{𝒃𝒆𝒂𝒎}\mathbf{}\mathbf{0.1}`$ by the $`𝐁`$ meson mass constraint, and it cannot exceed 1.0 due to the beam energy constraint.

After all these requirements, the distribution of the difference between the reconstructed energy and that of generated $`𝐁`$ hadrons in simulated data has a rms width of $`\mathbf{4.8}\mathrm{𝐆𝐞𝐕}`$. The energy dependence of the $`𝐁`$ hadron energy resolution is shown in Figure 2a. The complete $`𝐁`$ hadron selection applied to the full data sample results in $`\mathbf{270\hspace{0.17em}707}`$ tagged jets with a $`𝐛`$ purity of 96%. The average $`𝐁`$ hadron selection efficiency is 16%, with an energy dependence as shown in Figure 2b. The measured $`𝐁`$ hadron energy distribution, scaled to the beam energy, is shown in Figures 36, and compared to the various models described in the next section.

## 5 Test of hadronisation models

The $`𝐁`$ hadron energy distributions predicted by the JETSET 7.4, HERWIG 5.9, and HERWIG 6.2 Monte Carlo models are compared to the OPAL data. All Monte Carlo simulated events are passed through a detailed detector simulation . The comparison is performed using the distribution of the reconstructed scaled energy of the weakly decaying $`𝐁`$ hadron $`𝒙_𝐄`$.

The HERWIG Monte Carlo uses a parton shower fragmentation followed by cluster hadronisation model with few parameters. No parameters are varied in this analysis. This simplifies the model test to a mere comparison of the $`𝒙_𝐄`$ distributions obtained with data and Monte Carlo simulation. Both HERWIG versions are set up to conserve the initial $`𝐛`$ quark direction in the $`𝐁`$ hadron creation during cluster decay (cldir=1). The main difference between the two HERWIG samples used in this analysis is that Gaussian smearing of the $`𝐁`$ hadron direction around the initial $`𝐛`$ quark flight direction is applied in the HERWIG 5.9 sample (clsmr=0.35), while smearing is not used in the HERWIG 6.2 sample (clsmr(2)=0).

The JETSET Monte Carlo is based on a parton shower fragmentation followed by string hadronisation scheme. It requires a fragmentation function to describe the distribution of the fraction $`𝒛`$ of the string light cone momentum that is assigned to a hadron produced at the end of the string. The JETSET sample in this analysis is reweighted to use the fragmentation functions of Kartvelishvili et al. , Bowler , the Lund symmetric model , and the fragmentation functions of Peterson et al. , and Collins-Spiller . The Lund symmetric and Bowler functions are simplified by assuming the transverse mass of the $`𝐛`$ quark, $`𝒎_{\mathbf{}}`$, to be constant, which is justified by the smallness of the average transverse momentum compared to the $`𝐛`$ quark mass. A further simplification in the Bowler parametrisation is the assumption of an equality of $`𝐛`$ quark and hadron masses. The functional forms of the fragmentation functions are given in Table 2. The parameters of the respective fragmentation functions are fitted to obtain a best match of the observed $`𝒙_𝐄`$ distributions in data and Monte Carlo simulation. In the case of the Peterson et al., Collins-Spiller, and Kartvelishvili et al. models, one free parameter is available. The Lund and Bowler models each have two free fit parameters. A $`𝝌^\mathrm{𝟐}`$ fit is performed in 46 bins in the $`𝒙_𝐄`$ range of 0.5 to 0.95, where in all samples the number of candidates in each bin is large enough to justify the assumption of Gaussian errors on the bin content. The fragmentation function and its parameters are adjusted during the fit by reweighting the Monte Carlo simulated events, similar to the procedure applied in .

The reweighting fit is performed separately for each data sample, and the fit results are averaged with weights according to the size of the respective datasets. Consistent results are obtained for all datasets. The average parameter values are summarised in Table 3. For each parametrisation the corresponding model-dependent mean scaled energy of weakly decaying $`𝐁`$ hadrons is given. Data samples at $`\sqrt{𝒔}\mathbf{=}𝒎_𝐙\mathbf{+}\mathbf{1.8}\mathrm{𝐆𝐞𝐕}`$ are excluded in the calculation of the average $`\mathbf{}𝒙_𝐄\mathbf{}`$. Table 3 also gives a comparison of the fit quality of all JETSET 7.4 fits on the 1994 data, and the HERWIG 5.9 and HERWIG 6.2 results. The fit results on the 1994 data are shown in Figures 36. The ordering of the models according to the goodness of the fits in 1994 data agrees with all other large data samples; only in a few smaller samples is a slightly different ordering observed. The quoted $`𝝌^\mathrm{𝟐}`$/d.o.f. values only take into account the statistical uncertainty of data and Monte Carlo simulation. Systematic uncertainties are discussed later. The Bowler, Lund symmetric, and Kartvelishvili et al. models are preferred by the data, with respective $`𝝌^\mathrm{𝟐}`$/d.o.f. values of 67/44, 75/44, and 99/45 in the 1994 sample. Figures 36 show that the Peterson et al. and Collins-Spiller parametrisations for JETSET, as well as the HERWIG 6.2 model, are too broad. The HERWIG 5.9 model is too soft.

## 6 Model-independent measurement of $`\mathbf{}𝒙_𝐄\mathbf{}`$

In the previous section, information was extracted from the observed energy distribution making explicit use of a set of models to describe the data. In this section, a measurement of the mean scaled energy of $`𝐁`$ hadrons, $`\mathbf{}𝒙_𝐄\mathbf{}`$, outside a specific model framework will be presented. This is accomplished by unfolding the observed energy distribution.

Two complementary unfolding procedures are used to obtain an estimate of the true $`𝒙_𝐄`$ distribution from the observed distribution of the reconstructed scaled $`𝐁`$ hadron energy. In both cases the amount and energy distribution of background in the $`𝐁`$ hadron candidate sample is estimated from the Monte Carlo simulation and subtracted from the data.

The main method starts by fitting the observed data $`𝒙_𝐄`$ distribution, and the observed and the true $`𝒙_𝐄`$ distribution in the Monte Carlo simulation with smooth functions (splines). The true and observed Monte Carlo distributions are then reweighted simultaneously until the observed $`𝒙_𝐄`$ distribution agrees in data and simulation. The reweighted true $`𝒙_𝐄`$ distribution of the Monte Carlo simulation then provides an estimate of the corresponding distribution in data. Details of how the result is stabilised are described later. This method is almost independent of the initial Monte Carlo distribution and thus reduces model-dependence in the unfolding process. Furthermore, the result is represented as an unbinned spline function, which is optimal for the calculation of the mean value of the unfolded distribution. This algorithm is coded using the software package RUN  and was already used in .

The second approach makes use of the SVD-GURU software package . The correspondence between the observed and true $`𝐁`$ hadron energy distributions in the Monte Carlo simulation is represented by a $`\mathrm{𝟐𝟎}\mathbf{\times }\mathrm{𝟐𝟎}`$ matrix. The unfolding process comprises a matrix inversion to obtain an estimate of the true data $`𝒙_𝐄`$ distribution from the observed distribution. In this approach, the model dependence was found to be stronger than when using the RUN program. Furthermore, a coarse binning appropriately adapted to the detector’s resolution and the amount of available data might lead to systematic effects when describing the energy distribution in terms of its mean value. Therefore SVD-GURU is only used to cross-check the result obtained by RUN and to provide an estimate of the systematic uncertainty due to unfolding.

Raw unfolding solutions often oscillate strongly around the correct solution. In the case of a binned representation of the data this effect can simply be understood by strong negative bin-by-bin correlations introduced by the finite detector resolution. Both methods used here suppress these oscillations by limiting the number of degrees of freedom of the unfolding solution. The RUN algorithm represents the unfolding result as expansion into a set of orthogonal functions. The uncertainties on the coefficients of these functions are determined, and only those functions with coefficients significantly different from zero are taken into account. SVD-GURU rotates the unfolding matrix to estimate its effective rank. The unfolding is then performed in a rotated space with a smaller matrix including only the significant contributions. The number of degrees of freedom used for the unfolding procedure was found to agree in the RUN and SVD-GURU approaches in all subsamples described below. A further means of regularisation is available in the RUN package. Of all remaining solutions to the unfolding problem, one is chosen that minimises the integral over the squared first derivative of the unfolding solution. Monte Carlo studies show that this regularisation leads to essentially bias-free results on all samples. The performance of the unfolding algorithms is illustrated in Figure 7.

The unfolding is performed separately for data from all years of 1992–2000. The 1993 and 1995 datasets at a centre-of-mass energy below $`𝒎_𝐙`$ show $`𝒙_𝐄`$ distributions that are compatible with those taken at the $`𝐙`$ resonance peak in the Monte Carlo simulation. The 1993 and 1995 datasets at $`𝒎_𝐙\mathbf{+}\mathbf{1.8}\mathrm{𝐆𝐞𝐕}`$ show a significantly lower $`\mathbf{}𝒙_𝐄\mathbf{}`$, caused by a large amount of initial state radiation at this energy. In this case the quark energy prior to fragmentation is lower on average than the beam energy. As the beam energy is used as an estimator of the quark energy prior to fragmentation, the average $`𝒙_𝐄`$ value is lower than in samples without significant initial state radiation. The $`𝒎_𝐙\mathbf{+}\mathbf{1.8}\mathrm{𝐆𝐞𝐕}`$ samples are therefore analysed separately.

Both RUN and SVD-GURU analyses are performed with a Monte Carlo simulation that is reweighted to match the best result of the model-dependent reweighting fits for the respective datasets. This procedure is also followed by SLD in their latest $`𝐛`$ hadronisation analysis . The goal is to reduce the dependence of the unfolding result on the Monte Carlo sample used for unfolding. The effect of not using the best parametrisation, but the second and third best instead, are studied below as a systematic effect. JETSET 7.4 Monte Carlo simulation samples are used to obtain the central result, and will be compared to unfolding results using HERWIG in the discussion of systematic effects below.

The unfolding result for all data samples at or below the $`𝐙`$ resonance is shown in Figure 8. The mean scaled $`𝐁`$ hadron energy, obtained with the RUN unfolding algorithm and averaged over all datasets by using the subsample size for the weights, is

$$\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{0.7193}\mathbf{\pm }\mathbf{0.0016}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}\mathbf{,}$$

where the uncertainty includes the statistical uncertainties due to limited data and Monte Carlo sample sizes, and the statistical uncertainty on the Monte Carlo efficiency. Consistent $`\mathbf{}𝒙_𝐄\mathbf{}`$ values were obtained for the individual data samples (see Figure 8). This model-independent measurement agrees with the $`\mathbf{}𝒙_𝐄\mathbf{}`$ values in the framework of the best models as seen in Table 3. The unfolding result spline for the full dataset is plotted in Figure 9.

The mean scaled energy observed in the $`𝒎_𝐙\mathbf{+}\mathbf{1.8}\mathrm{𝐆𝐞𝐕}`$ samples is found to be

$$\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{0.7130}\mathbf{\pm }\mathbf{0.0056}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}\mathbf{.}$$

The 1993 and 1995 $`𝒎_𝐙\mathbf{+}\mathbf{1.8}\mathrm{𝐆𝐞𝐕}`$ data samples give consistent values of $`\mathbf{}𝒙_𝐄\mathbf{}`$. The difference of the results for the different energies is consistent with the prediction obtained from Monte Carlo samples at similar energies.

The results obtained with the unfolding program SVD-GURU ($`\mathbf{0.7195}\mathbf{\pm }\mathbf{0.0015}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}`$ for the main dataset, $`\mathbf{0.7152}\mathbf{\pm }\mathbf{0.0053}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}`$ for the $`𝒎_𝐙\mathbf{+}\mathbf{1.8}\mathrm{𝐆𝐞𝐕}`$ samples) are in very good agreement with the ones achieved with the RUN package. This is also demonstrated in Figure 9, where the results obtained from the full $`\sqrt{𝒔}\mathbf{}𝒎_𝐙`$ data sample by both algorithms are compared. The statistical uncertainties of both methods are also very similar.

## 7 Systematic uncertainties

Given the large data sample collected with the OPAL detector, and the inclusive character of the analysis presented here, the statistical uncertainties on the results of the previous sections are expected to be small compared to the systematic uncertainties introduced by limited knowledge of physics parameters which possibly affect the measured quantities. In this section, an overview of all systematic checks is given for the reweighting fit and for the unfolding analysis.

The distribution of the $`𝐛`$-tagging discriminant in data agrees well with the Monte Carlo simulation. However, it is necessary for this analysis to ensure that this is separately true for different $`𝐁`$ hadron energy regions. Therefore the $`𝐛`$-tagging discriminant was investigated in ten bins in $`𝒙_𝐄`$. The ratio of the $`𝐛`$-tagging discriminant distributions of data and Monte Carlo simulation is calculated and fitted by a linear function in each energy bin separately. The slope of this function is an indicator of the quality of agreement of data and simulation. All fitted slopes are compatible with zero within two standard deviations.

### Systematic uncertainties for fits to various models

The systematic checks performed for the fits to models are listed below. The resulting systematic uncertainty estimates for the mean scaled energy $`\mathbf{}𝒙_𝐄\mathbf{}`$ in the framework of the respective models are summarised in Table 4.

* The energy resolution of the OPAL calorimeters in the Monte Carlo simulation is varied by $`\mathbf{\pm }\mathrm{𝟏𝟎}\mathbf{\%}`$ relative to its central value. This range is motivated by jet energy resolution studies in two-jet events, where a difference between the resolution in data and the Monte Carlo simulation of 3.6% was found in some datasets. Under the assumptions that 50% of the total jet energy are contributed by neutral particles and that the observed difference can be fully accounted to the calorimeters, the above variation range covers this effect.
* Similar studies indicate a possible difference of the energy scale between data and simulation of up to 0.4% in some datasets. The energy scale is varied within this range, and the resulting difference of the fit results is taken as systematic uncertainty.
* The resolution of the $`𝒓\mathit{\varphi }`$-related track parameters $`𝒅_\mathrm{𝟎}`$ (transverse impact parameter), $`\mathit{\varphi }_\mathrm{𝟎}`$ (initial azimuth), and $`𝜿`$ (curvature) is varied in the range $`\mathbf{\pm }\mathrm{𝟏𝟎}\mathbf{\%}`$ .
* The resolution of the $`𝒓𝒛`$-related track parameters $`𝒛_\mathrm{𝟎}`$ (longitudinal impact parameter) and $`\mathrm{𝐭𝐚𝐧}𝝀`$ (dip angle, $`𝝀\mathbf{=}𝝅\mathbf{/}\mathrm{𝟐}\mathbf{}𝜽`$) is varied by $`\mathbf{\pm }\mathrm{𝟏𝟎}\mathbf{\%}`$ .
* Figure 2a shows a large energy reconstruction bias for low energy $`𝐁`$ hadrons. Both the number of low energy $`𝐁`$ hadrons and the efficiency for reconstructing them are small. This large bias therefore only affects a very small fraction of the candidates. The effect of possible mismodeling of the bias in the Monte Carlo simulation is evaluated by varying the bias around the values that describe the data best. The reconstruction bias in the high energy region cannot be changed by more than $`\mathbf{\pm }\mathrm{𝟏}\mathbf{\%}`$ without leading to a significant degradation of the agreement between data and Monte Carlo simulation. The bias for low energy $`𝐁`$ hadrons, with a reconstructed $`𝒙_𝐄`$ below 0.6, is varied by $`\mathbf{\pm }\mathrm{𝟏𝟎}\mathbf{\%}`$. This range is also motivated by significant degradation of the agreement between data and Monte Carlo simulation. The largest deviations of the measured quantities are taken as systematic uncertainties.
* The relative fraction of different $`𝐁`$ hadron species in the sample of primary $`𝐁`$ hadrons influences the measurement, because different $`𝐁`$ hadron species have different energy distributions. All values obtained in this analysis are calculated with Monte Carlo samples that are reweighted to reflect the current best knowledge of the hadron fractions. The associated systematic uncertainty is estimated by varying the $`𝐛`$ baryon fraction within the range $`\mathbf{(}\mathbf{10.3}\mathbf{\pm }\mathbf{1.8}\mathbf{)}\mathbf{\%}`$  given by the average of the LEP/SLD/CDF measurements of this quantity. The fraction of $`𝐁_𝐬`$ in the sample is varied in the range of $`\mathbf{(}\mathbf{9.8}\mathbf{\pm }\mathbf{1.2}\mathbf{)}\mathbf{\%}`$ .
* The amount of orbitally excited $`𝐁_𝑱^{\mathbf{(}\mathbf{}\mathbf{)}}`$ mesons has been measured by all LEP collaborations . An error weighted average of the LEP measurements is $`\mathbf{(}\mathbf{28.4}\mathbf{\pm }\mathbf{3.5}\mathbf{)}\mathbf{\%}`$, and the fraction of orbitally excited $`𝐁_𝑱^{\mathbf{(}\mathbf{}\mathbf{)}}`$ mesons is varied within this range.
* The Q-value of orbitally excited $`𝐁_𝑱^{\mathbf{(}\mathbf{}\mathbf{)}}`$ mesons  is about $`\mathrm{𝟒𝟎}\mathrm{𝐌𝐞𝐕}`$ smaller than the value in the Monte Carlo samples used in this analysis. All results are corrected for this effect, and the difference to the values obtained without correction is taken as systematic uncertainty.
* The average multiplicity of charged particles from a $`𝐁`$ hadron decay was found at LEP to be $`\mathbf{4.955}\mathbf{\pm }\mathbf{0.062}`$ , and is varied within this range.
* The average lifetime of weakly decaying $`𝐁`$ hadrons affects the efficiency of secondary vertex reconstruction, and is varied in the range $`\mathbf{(}\mathbf{1.577}\mathbf{\pm }\mathbf{0.016}\mathbf{)}\mathrm{𝐩𝐬}`$ .
* The average lifetime of weakly decaying charm hadrons determines the amount of charm background found in the $`𝐁`$ hadron candidate sample. The $`𝐃^\mathrm{𝟎}`$, $`𝐃^\mathbf{+}`$, $`𝐃_{𝐬}^{}{}_{}{}^{\mathbf{+}}`$, and $`𝚲_{𝐜}^{}{}_{}{}^{\mathbf{+}}`$ lifetimes are varied within the uncertainties quoted in .
* The charm background in the Monte Carlo simulation samples is reweighted to the same hadronisation model as the $`𝐁`$ hadron distribution in the respective fits. The central values and uncertainties of the charm fragmentation function are taken from earlier OPAL measurements . For the evaluation of the systematic uncertainty the parameters are varied within their uncertainties.
* Jets from gluon splitting to $`𝐛\overline{𝐛}`$ quark pairs are treated as background in this analysis. To account for the uncertainty of the average LEP result of the gluon splitting rate, $`\mathbf{0.00254}\mathbf{\pm }\mathbf{0.00050}`$ $`𝐛\overline{𝐛}`$ pairs per hadronic event , the rate is varied within this range.
* The number of $`𝐜\overline{𝐜}`$ pairs from gluon splitting per hadronic event is varied within the LEP uncertainty of $`\mathbf{0.0299}\mathbf{\pm }\mathbf{0.0039}`$ .
* The partial width of the $`𝐙`$ into $`𝐛\overline{𝐛}`$ quark pairs, normalised to the total hadronic width of the $`𝐙`$, is measured to be $`𝐑_𝐛\mathbf{=}\mathbf{0.21646}\mathbf{\pm }\mathbf{0.00065}`$ . Varying this fraction within the quoted uncertainty leads to varying background levels in the unfolding Monte Carlo sample. This causes negligible changes of the fit results.
* The analogous quantity for charm quark pairs, $`𝐑_𝐜`$, is less well known, with a current best value of $`\mathbf{0.1719}\mathbf{\pm }\mathbf{0.0031}`$ . The impact on the fit results of varying $`𝐑_𝐜`$ within this range is negligible.
* Limited knowledge of the LEP beam energy produces a corresponding uncertainty on $`𝒙_𝐄`$, although the dependency of $`𝒙_𝐄`$ on the beam energy is reduced due to the fact that the beam energy also enters the calculation of the reconstructed $`𝐁`$ hadron energy via the beam energy constraint. The assumed LEP beam energy is varied within $`\mathbf{\pm }\mathrm{𝟖}\mathrm{𝐌𝐞𝐕}`$, which is the largest reported uncertainty for any sample at or close to the $`𝐙`$ resonance . A correlation of 100% between the resulting uncertainties for the different data taking years is assumed.
* The parameter values depend slightly on the $`𝒙_𝐄`$ range used for the fit. The lower end of the fit range is varied within $`𝒙_𝐄\mathbf{=}\mathbf{0.5}\mathbf{\pm }\mathbf{0.1}`$, and the upper range is varied within $`𝒙_𝐄\mathbf{=}\mathbf{0.95}\mathbf{\pm }\mathbf{0.05}`$. The largest deviation from the result obtained using the central values is taken as the systematic uncertainty.
* The bin width used in the fit is varied by $`\mathbf{\pm }\mathrm{𝟏𝟎}\mathbf{\%}`$. The maximum deviation from the result with standard binning is used to estimate the associated systematic uncertainty.

The JETSET 7.4 Monte Carlo samples used for the reweighting fit were generated using the Peterson et al. fragmentation function with $`𝜺_𝐛\mathbf{=}\mathrm{𝟑𝟖}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟒}`$. The fit result for the Peterson et al. function in data is $`𝜺_𝐛\mathbf{=}\mathbf{(}\mathbf{41.2}\mathbf{\pm }\mathbf{0.7}\mathbf{)}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟒}`$. The fact that the Monte Carlo tuning and the data fit result are close has the advantage that adverse effects due to weights far from 1.0 are not expected. However, additional studies were performed to verify that the closeness of the two values is not introduced by improper reweighting. An older sample of 4 million hadronic JETSET 7.4 events with Peterson et al. fragmentation function with $`𝜺_𝐛\mathbf{=}\mathrm{𝟓𝟕}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟒}`$ is used to repeat the fit for the 1994 data sample. The fit result obtained with this sample ($`𝜺_𝐛\mathbf{=}\mathbf{(}\mathbf{40.3}\mathbf{\pm }\mathbf{1.1}\mathbf{)}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟒}`$) is in agreement with the result obtained with the main $`𝜺_𝐛\mathbf{=}\mathrm{𝟑𝟖}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟒}`$ 1994 Monte Carlo sample ($`𝜺_𝐛\mathbf{=}\mathbf{(}\mathbf{40.6}\mathbf{\pm }\mathbf{1.0}\mathbf{)}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟒}`$).

### Systematic uncertainties of the unfolding analysis

The same systematic uncertainties as for the fragmentation function fits are evaluated also for the $`\mathbf{}𝒙_𝐄\mathbf{}`$ measurement, with the exception of fit range effects, which are specific to the model-dependent fit procedure. Binning effects are not present in RUN. The resulting systematic uncertainties are summarised in Table 5.

An additional uncertainty is introduced by the dependence of the $`\mathbf{}𝒙_𝐄\mathbf{}`$ measurement on detector and acceptance modelling in the Monte Carlo simulation. As in all unfolding problems, one needs the resolution (or spectral) function $`𝒈\mathbf{(}𝒙_\mathrm{𝟎}\mathbf{,}𝒙\mathbf{)}`$ where $`𝒙_\mathrm{𝟎}`$ is the energy of the hadrons entering the detector and $`𝒙`$ their measured energy. This function is not measured, but calculated by the OPAL detector simulation . Since the detector simulation is generally made in the framework of some specific Monte-Carlo program generating hadron distributions, a small residual dependency of $`𝒈\mathbf{(}𝒙_\mathrm{𝟎}\mathbf{,}𝒙\mathbf{)}`$ on the particular Monte Carlo event generator remains. To estimate the associated systematic uncertainty, the unfolding procedure was repeated using not only the best reweighting fit result, but also all other models as initial estimators of the true distribution. This study was performed independently for all datasets. The unfolding was also performed with the JETSET Monte Carlo sample replaced by HERWIG 5.9 and 6.2 samples. The check using JETSET with the Collins-Spiller parametrisation dominates the modelling uncertainty in the negative direction, which is taken as the largest observed deviation from the central $`\mathbf{}𝒙_𝐄\mathbf{}`$ value. The uncertainty in the positive direction is dominated by the third best model, which in almost all datasets is the Kartvelishvili et al. parametrisation. The resulting model uncertainty is found to be $`{}_{\mathbf{}\mathbf{0.0016}}{}^{}{}_{}{}^{\mathbf{+}\mathbf{0.0024}}`$.

The result of the unfolding procedure with the RUN algorithm is cross-checked with the SVD-GURU method, and the difference between the two results is assigned as systematic uncertainty. Furthermore, a difference of similar size is observed between the spline unfolding result of the RUN method and a binned representation of the unfolded distribution. This difference is also included in the unfolding method uncertainty.

Summing all systematic uncertainties in quadrature yields a total systematic uncertainty on $`\mathbf{}𝒙_𝐄\mathbf{}`$ of $`{}_{\mathbf{}\mathbf{0.0033}}{}^{}{}_{}{}^{\mathbf{+}\mathbf{0.0038}}`$. As expected, the systematic uncertainty is larger than the statistical precision.

Table 6 shows a representation of the RUN result in 20 bins in the $`𝒙_𝐄`$ range between 0.1 and 1.0. This table, along with the full correlation matrix of statistical (Table 7) and systematic (Tables 8 and 9) uncertainties can be used to compare the OPAL results with further hadronisation models not discussed here. It has to be pointed out again that the $`\mathbf{}𝒙_𝐄\mathbf{}`$ value obtained from the binned RUN result is smaller than the unbinned result by $`𝚫\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{}\mathbf{0.0002}`$. This is caused by binning effects which are reduced by using a small bin width, but cannot be entirely avoided.

## 8 Summary and discussion

Using an unfolding technique to reduce the dependence on $`𝐛`$ quark hadronisation models, the mean scaled energy of weakly decaying $`𝐁`$ hadrons in $`𝐙`$ decays is measured to be

$$\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{0.7193}\mathbf{\pm }\mathbf{0.0016}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}_{\mathbf{}\mathbf{0.0033}}^{\mathbf{+}\mathbf{0.0038}}\mathbf{(}𝒔𝒚𝒔𝒕\mathbf{)}\mathbf{.}$$

This is the most precise available measurement of this quantity. Consistent results are obtained using an alternative unfolding method and from model-dependent reweighting fits.

The result obtained here is in good agreement with a recent result from the ALEPH Collaboration , $`\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{0.716}\mathbf{\pm }\mathbf{0.006}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}\mathbf{\pm }\mathbf{0.006}\mathbf{(}𝒔𝒚𝒔𝒕\mathbf{)}`$. ALEPH uses exclusive semileptonic $`𝐁`$ decays, leading to a smaller candidate sample and thus a larger statistical uncertainty. Another new measurement by SLD  gives a somewhat lower value: $`\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{0.709}\mathbf{\pm }\mathbf{0.003}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}\mathbf{\pm }\mathbf{0.003}\mathbf{(}𝒔𝒚𝒔𝒕\mathbf{)}\mathbf{\pm }\mathbf{0.002}\mathbf{(}𝒎𝒐𝒅𝒆𝒍\mathbf{)}`$. The difference between the OPAL and the SLD measurements has a statistical significance of less than $`\mathrm{𝟐}`$ standard deviations taking only the uncorrelated uncertainties into account. Another $`\mathbf{}𝒙_𝐄\mathbf{}`$ measurement was recently performed using inclusive $`𝐁\mathbf{}\mathbf{}`$ decays . Modelling of the lepton energy spectrum introduces additional systematic errors in the lepton-based analysis. The result of $`\mathbf{}𝒙_𝐄\mathbf{}\mathbf{=}\mathbf{0.709}\mathbf{\pm }\mathbf{0.003}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}\mathbf{\pm }\mathbf{0.003}\mathbf{(}𝒔𝒚𝒔𝒕\mathbf{)}\mathbf{\pm }\mathbf{0.013}\mathbf{(}𝒎𝒐𝒅𝒆𝒍\mathbf{)}`$ is compatible with the analysis presented here, especially given that the result in  is not model-independent, but based on a Peterson et al. parametrisation. The LEP average result for $`\mathbf{}𝒙_𝐄\mathbf{}`$ in the framework of the Peterson et al. model, obtained from earlier analyses , is $`\mathbf{0.702}\mathbf{\pm }\mathbf{0.008}`$, again in excellent agreement with the value of $`\mathbf{0.7023}\mathbf{\pm }\mathbf{0.0006}\mathbf{(}𝒔𝒕𝒂𝒕\mathbf{)}\mathbf{\pm }\mathbf{0.0019}\mathbf{(}𝒔𝒚𝒔𝒕\mathbf{)}`$ found in this analysis.

The best description of the data with a fragmentation function with one free parameter is achieved with the Kartvelishvili et al. model. The Peterson et al. and Collins-Spiller models produce energy distributions which are too broad to describe the data. Similar features have been observed by SLD and ALEPH in their recent publications. The Bowler and Lund parametrisations, each having two free parameters, achieve a better $`𝝌^\mathrm{𝟐}`$/d.o.f. in this analysis and are clearly compatible with the data. The same conclusion is reached by SLD while ALEPH did not test these models. The HERWIG cluster model is clearly disfavoured. The main difference of the two HERWIG versions tested in this analysis is the amount of smearing of the $`𝐁`$ hadron direction around the initial $`𝐛`$ quark direction. Significant smearing is used in the HERWIG 5.9 sample, softening the spectrum too much. The HERWIG 6.2 sample is used without any smearing, giving an $`𝒙_𝐄`$ distribution which is in much better agreement with the data, but which is still too broad. Similar results are obtained by SLD.

The fitted values of the parameters describing each hadronisation model agree less well between the different experiments than the measured $`\mathbf{}𝒙_𝐄\mathbf{}`$ values. The parameter values depend critically on details of the Monte Carlo tuning, which is not identical in all respects among the collaborations, although efforts have been made to correct most relevant Monte Carlo parameters to a common set of values.

A general conclusion of the analysis presented here is that the parton shower plus string hadronisation Monte Carlo models provide a good description of the current data. The fragmentation functions derived from intrinsic symmetries of the string model (Bowler, Lund symmetric) are favoured over the phenomenological approaches of Kartvelishvili et al., Peterson et al., and Collins-Spiller.

Acknowledgements

We particularly wish to thank the SL Division for the efficient operation of the LEP accelerator at all energies and for their close cooperation with our experimental group. In addition to the support staff at our own institutions we are pleased to acknowledge the
Department of Energy, USA,
National Science Foundation, USA,
Particle Physics and Astronomy Research Council, UK,
Natural Sciences and Engineering Research Council, Canada,
Israel Science Foundation, administered by the Israel Academy of Science and Humanities,
Benoziyo Center for High Energy Physics,
Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) and a grant under the MEXT International Science Research Program,
Japanese Society for the Promotion of Science (JSPS),
German Israeli Bi-national Science Foundation (GIF),
Bundesministerium für Bildung und Forschung, Germany,
National Research Council of Canada,
Hungarian Foundation for Scientific Research, OTKA T-029328, and T-038240,
Fund for Scientific Research, Flanders, F.W.O.-Vlaanderen, Belgium.


format: text

text:
# Angular momenta creation in relativistic electron-positron plasma

## I Introduction

The problem of electromagnetic (EM) wave propagation and related phenomena in relativistic plasmas has attracted considerable attention in the recent past. From the nonthermal emission of the high-energy radiation coming from a variety of compact astrophysical objects it has become possible to deduce the presence of a population of relativistic electrons in the plasma created in the dense radiation fields of those sources 1 . The principal components of these plasmas could be either relativistic electrons and nonrelativistic ions (protons), or relativistic electron-positron (e-p) pairs.

Relativistic e-p dominated plasmas may be created in a variety of astrophysical situations. The e-p plasmas are likely to be found in pulsar magnetospheres pulsar , in the bipolar outflows (jets) in Active Galactic Nuclei (AGN) AGN , and at the center of our own Galaxy Galaxy . The presence of e-p plasma is also argued in the MeV epoch of the early Universe. In the standard cosmological model, temperatures in the MeV range ($`T10^{10}`$K$`1`$MeV) prevail up to times $`t=1`$sec after the Big Bang universe . In this epoch, the main constituent of the Universe is the relativistic e-p plasma in equilibrium with photons, neutrinos, and antineutrinos.

Contemporary progress in the development of super strong laser pulses with intensities $`I10^{2123}\mathrm{W}/\mathrm{cm}^2`$ has also made it possible to create relativistic plasmas in the laboratory by a host of experimental techniques laser1 . At the focus of an ultrastrong laser pulses, the electrons can acquire velocities close to the speed of light opening the possibility of simulating in the laboratory the conditions and phenomena that, generally, belong in the astrophysical realm laser2 .

Elucidation of the electromagnetic wave dynamics in a relativistic plasmas will, perhaps, be an essential tool for understanding the radiation properties of astrophysical objects as well as of the media exposed to the field of superstrong laser radiation. Although the study of wave propagation in relativistic plasmas has been in vogue for some time, it is only in the recent years that the nonlinear dynamics of EM radiation in e-p dominated plasmas 8  has come into focus. The enhanced interest stems from two facts: 1) e-p plasmas seem to be essential constituents of the universe, and 2) under certain conditions, even an ultrarelativistic electron-proton plasma can behave akin to an e-p plasma 9 .

Recently, we firstly found dark soliton as well as vortex soliton solutions in e-p plasma TBM-01 . In Ref. FaBu-dark , it is also shown analytically that dark soliton is the natural nonlinear coherent structure in unmagnetized cold e-p plasmas. However, it is conceivable that soliton solutions obtained in a one dimensional formulation will turn out to be unstable in higher dimensions, which may lead to the creation of vortex solitons. Since dark and vortex solitons are asymptotically nonvanishing, they have received much less attention than their localized cousins due to the generally accepted requirement that the fields be localized in physical system. However, in recent experiments studying laser field dynamics in different kinds of optical media, it was demonstrated that dark and vortex solitons can be readily created as superimpositions upon a localized field background KL-98 . Vortex soliton solutions are also found in imperfect Bose gas in the superfluids Pitaevskii , and they are extensively investigated and discussed BEcondense . In this paper, we systematically investigate the instability of dark solitons and show that it can lead to the formation of vortex solitons. We speculate about interesting application of the vortex soliton in the early universe.

In the recent paper TBM-01 , we had developed an argument for the creation of domains of nonzero angular momentum in the MeV era of the early universe when it is supposed to be dominated by a plasma of e-p pairs. We first showed that in such a plasma the dynamics of a pulse of electromagnetic radiation, with a frequency much larger than the plasma frequency, is controlled by a generalized nonlinear Schrödinger equation (GNSE) with a defocusing nonlinearity. Then, borrowing a result of nonlinear optics KL-98  where the standard nonlinear Schrödinger equation (NSE) with a cubic nonlinearity has been investigated in great depth and detail, we conjectured that even the GNSE whose nonlinearity is similar in nature to that of NSE will allow vortex soliton solutions with an angular momentum that is conserved during propagation. The latter system also allows dark soliton solutions in one dimension which are known to be unstable to two-dimensional (2d) perturbations and eventually evolve to 2d vortex chains.

In this paper we demonstrate that the arguments given in Ref. TBM-01  are actually borne out by the direct solutions of the derived GNSE. We begin by analytically showing that the 1d dark solitons of the GNSE are, indeed, unstable to 2d perturbations. In the process, we derive the instability criterion which agrees with the numerical solution of the linearized system.

The main part of the demonstration, however, comes from a numerical simulation of our GNSE. Starting from a broad variety of initial conditions, we find the eventual emergence of angular momentum carrying vortex solitons. That is, the stipulated solitonic structures are readily accessible within the framework of this general equation with a nonlinearity more complicated than that of NSE. Thus we can state with much greater confidence that electromagnetism, operative in the MeV era, could easily be the primordial source of angular momentum associated with various structures of the observable universe.

We again stress that the present results are quite general in the nonlinear dynamics in an electron positron plasma.

## II Single vortex soliton

In the envelope approximation, a finite amplitude, circularly polarized electromagnetic (EM) pulse propagating in a relativistic e-p plasma obeys (for details, see TBM-01 ):

$$2\mathrm{i}\omega _\tau A_{}+c^2_{}^2A_{}+\frac{c^2}{\gamma _g^2}_\xi ^2A_{}+\frac{2\omega _e^2}{G_0}\left(1\frac{G_0}{\gamma G}\frac{n}{n_0}\right)A_{}=0,$$
(1)
where $`A_{}`$ is the slowly varying amplitude of the perpendicular (to the propagation direction) vector potential, $`_{}^2=_x^2+_y^2`$, and $`\xi =zv_gt`$ is the “comoving” (with group velocity $`v_g`$) coordinate. Here $`\omega _e=(4\pi n_{0e}e^2/m_e)^{1/2}`$ is the electron Langmuir frequency, $`n_0=n_{0e}=n_{0p}`$ is the number density of the unperturbed background (subscript $`0`$ denotes the value at infinity), $`\gamma `$ is the Lorentz factor, and $`G=K_3(m_ec^2/T)/K_2(m_ec^2/T)`$ with $`K_n`$, denoting the $`n`$-th order modified Bessel functions of the second kind, is the enthalpy density.

In this approximation the continuity equation, determining the number density $`n`$, becomes

$$\frac{1}{\gamma G}\frac{n}{n_0}=\frac{v_g\gamma _g/c}{[\gamma _g^2G_0^2G^2e^2|A_{}|^2/(m_ec^2)^2]^{1/2}},$$
(2)
and the system is closed by the adiabatic equation of state

$$\frac{nm_ec^2/n_0T}{\gamma K_2(m_ec^2/T)}\mathrm{exp}(m_ec^2G/T)=\text{const.}$$
(3)
The low frequency motion of the plasma is driven by the ponderomotive pressure \[$`(𝒑_{e,p})^2`$\] of the high-frequency EM field, and is independent of the sign of the particles’ charge. It is perfectly natural to assume that the electron and the positron fluids have equal temperatures ($`T_{0e,p}=T_0`$) in equilibrium so that their effective masses ($`G_{e,p}=G`$) will also be equal. The radiation pressure will impart equal low-frequency momenta to both fluids allowing the possibility of overall density changes without producing charge separation. The charge neutrality conditions $`n_e=n_p=n`$, $`\varphi =0`$ has been assumed by neglecting the small inequality of the charge due to baryon asymmetry. It is also evident that the symmetry between the two fluids keeps their temperatures always equal ($`T_{e,p}=T`$) if they were equal initially. In deriving (1)-(2), we have also assumed that the plasma is transparent (i.e. $`\omega \omega _e`$), and that the longitudinal extent of the pulse is much shorter than its transverse dimensions ($`L_{}L_{}`$).

Defining the normalized variables

$$\frac{\omega _e^2}{2\omega G_0}tt,\frac{\omega _e}{c\sqrt{2G_0}}𝒓𝒓,\frac{e}{m_ec^2}\frac{1}{\gamma _gG_0}A_{}A_{},$$
(4)
and noting that even with the assumption $`_\xi _{}`$, the diffractive term can be the same order or even greater than the dispersive one for a highly transparent plasma ($`\gamma _g1`$), Eq. (1) converts to the following GNSE ($`\gamma _gG_0G`$):

$$\mathrm{i}_tA+\frac{1}{2}_{}^2A+f(|A|^2)A=0,$$
(5)
with the nonlinearity

$$f(|A|^2)=2\left(\frac{1}{\sqrt{1|A|^2}}1\right).$$
(6)
Equation (5) admits a symmetric two-dimensional solitary wave solution. For stationary solitons, the ansatz

$$A=A(r)\mathrm{exp}(\mathrm{i}m\theta \mathrm{i}\lambda t),$$
(7)
reduces (5) to

$$\frac{\mathrm{d}^2}{\mathrm{d}r^2}A+V^{}(A)=\frac{1}{r}\frac{\mathrm{d}A}{\mathrm{d}r}+\frac{m^2}{r^2}A,$$
(8)
with the potential $`V(A)`$ defined as

$$V(A)=(\lambda +2)A^2+4\sqrt{1A^2}4,$$
(9)
with the prime denoting the derivative with respect to its argument, and $`\lambda `$ representing the nonlinear frequency shift.

A numerical analysis of the solution of this eigenvalue problem is given in Ref. TBM-01 . Because of the absence of $`\xi `$-derivatives in (5), any pulse-like localized structure is allowed in the propagation direction. In the direction perpendicular to the propagation, on the other hand, the amplitude was shown to approach a constant value. Since there does exist a small baryon asymmetry in the early Universe that acts on a longer characteristic length BM-95 ; MBM-98 , the extent of the constant amplitude region may be considered to be finite.

Due to the single-valuedness of the vector potential, $`m`$ must be an integer and $`A`$ must vanish at the origin for non-zero values of $`m`$. The nonzero $`m`$ solutions are particularly important because they carry the orbital angular momentum $`𝑴`$

$$(𝑴)_z=\frac{\mathrm{i}}{2}\mathrm{d}𝒓_{}[𝒓_{}\times (A^{}_{}A\mathrm{c}.\mathrm{c}.)]_z.$$
(10)
It is straightforward to show that the system (5) conserves angular momentum, and the expression (10) is just the paraxial approximation for the orbital angular momentum, $`𝑴_E=d𝒓[𝒓\times (𝑬\times 𝑩)]`$, of the EM field AB-92 . The angular momentum carried by the vortex is $`M_z=mN`$ where the photon number $`N=d𝒓_{}|A|^2`$ is another conserved quantity; $`m`$ is also known as the “topological charge.” Strictly speaking, one must redefine the integrals of motion for non-vanishing boundary conditions perturbation theory , but such a renormalization is not important here because of the fact that infinite-extent solution is just a formal approximation. The presence of a small fraction of ions makes the physical solution decay at infinity BM-95 ; MBM-98 .

## III Dynamics of angular momenta creation

It is already suggested in Ref. TBM-01  that dark stripe solitons are unstable and break into vortex filaments. When the amplitude is small ($`|A|1`$), the nonlinearity reduces to the simple Kerr-type, and the vortex dynamics of NSE with such a nonlinearity has been studied in a variety of numerical calculations KL-98 ; perturbation theory ; McDonald ; Iv-97 . Our full nonlinearity is somewhat different and it has to be independently investigated. For this section, our goal is to show that the one-dimensional dark stripe solution is unstable against transverse perturbations in two dimensions. Such an instability causes the breakup of the stripe leading to a chain of vortex solitons with alternating polarity.

### III.1 Transverse instability

The one-dimensional stationary solution for the GNSE (5) has the form $`g(x)e^{\mathrm{i}\lambda t}`$ ($`\lambda >0`$), where the real function $`g(x)`$ satisfies

$$\lambda g+\frac{1}{2}\frac{\mathrm{d}^2g}{\mathrm{d}x^2}+f(g^2)g=0.$$
(11)
Since we are examining the linear stability of the solution with the least number of nodes, $`g`$ is assumed to be an odd function with a single node at the origin.

Let us perturb the one-dimensional solution by

$$A=(g+u+\mathrm{i}v)e^{\mathrm{i}\lambda t},$$
(12)
where the small perturbations $`u`$ and $`v`$ are real functions of $`x`$ and $`y`$.

Assuming sinusoidal behavior, $`u\mathrm{cos}(ky\mathrm{\Omega }t)`$ and $`v\mathrm{sin}(ky\mathrm{\Omega }t)`$, the linearized GNSE could be reduced to the following eigenvalue problem:

$$\mathrm{\Omega }^2v=_1_0v\frac{k^2}{2}(_0+_1)v+\frac{k^4}{4}v.$$
(13)
where

$`_0`$ $`=`$ $`\lambda +{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{x^2}}+f(g^2),`$ (14)
$`_1`$ $`=`$ $`\lambda +{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{x^2}}+f(g^2)+2g^2f^{}(g^2).`$ (15)
Remembering that $`_0g=0`$ and $`_1g_x=0`$, we can construct from (13) the relation

$$\mathrm{\Omega }^2g_x|v=\frac{k^4}{4}g_x|v+g_x|_1_0v\frac{k^2}{2}g_x|\mathrm{\hspace{0.17em}2}_1v\frac{k^2}{2}g_x|2g^2f^{}v$$
(16)
by multiplying $`g_x`$ on both sides and integrating with respect to $`x`$. Self-adjointness of $`_1`$ reduces (16) to

$$\mathrm{\Omega }^2=\frac{k^4}{4}\frac{k^2}{2}\frac{g_x|2g^2f^{}v}{g_x|v},$$
(17)
provided $`g_x|v`$ exists. A necessary condition for an exponential instability, then, is

$$\frac{g_x|2g^2f^{}v}{g_x|v}>0.$$
(18)
If it is true, then for sufficiently small $`|k|`$, we would have an instability, i.e., (18) is also a sufficient condition. Moreover, it is explicitly shown that $`\omega 0`$ in the limit $`k0`$ unless $`g_x|v`$ diverges.

We now make an approximate estimate of the ratio on the left hand side of (18). Let us assume that the operators $`_0`$ and $`_1`$ are of order unity and $`\mathrm{\Omega }^2`$ and $`k^2`$ are, in some sense, small. In this limit, we could estimate $`v`$ by solving

$$_1_0v=0$$
(19)
which has the solution $`v=g+q`$, where

$$_0q=g_x.$$
(20)
Notice that for $`v=g`$, the integrals occurring in (17) are zero and this part will not contribute to the integral; only $`q`$ will. From (20), we can deduce the following: the first consequence

$`g_x|g_x`$ $`=`$ $`g_x|_0q=g_x|(_0_1)q`$ (21)
$`=`$ $`g_x|2g^2f^{}q=g_x|2g^2f^{}v`$
converts (17) into

$$\mathrm{\Omega }^2=\frac{k^4}{4}\frac{k^2}{2}\frac{g_x|g_x}{g_x|v}.$$
(22)
To derive the second, we notice that

$$g_x|v=g_x|q=g_x|_0^1g_x$$
(23)
In principle, one can evaluate (23) in a rather straightforward way. But a very approximate estimate can be made by simply dropping $`_x^2`$ in $`_0`$ so that

$$g_x|v\frac{}{}g_x|\frac{g_x}{\lambda +f}$$
(24)
and

$$\mathrm{\Omega }^2=\frac{k^4}{4}\frac{k^2}{2}\frac{g_x|g_x}{g_x|g_x/(\lambda +f)}\frac{k^4}{4}\alpha k^2.$$
(25)
In (25) we have overestimated the denominator. The eigenvalue is

$$\mathrm{\Omega }=\mathrm{i}\alpha k\left[1\frac{k^2}{4\alpha }\right]^{1/2},$$
(26)
implying a window $`0<k<2\alpha ^{1/2}`$ in $`k`$, where instability is possible . From the property that the one-dimensional solution is stable to one-dimensional perturbation, the growth rate $`\mathrm{\Omega }`$ approaches zero as the wave number $`k`$ tends to zero.

In order to accomplish the numerical analysis, we introduce a complex-valued function $`w=u+\mathrm{i}v`$ in (12). Then, the linearized equation looks

$$\mathrm{i}_tw=\lambda w+\frac{1}{2}_{}^2w+f(g^2)w+2g^2f^{}(g^2)Re(w),$$
(27)
where the prime denotes derivative with respect to its argument. Putting $`_t=\mathrm{i}\mathrm{\Omega }`$ and $`_y=\mathrm{i}k`$ yields a linear eigenvalue equation for $`w`$:

$$\mathrm{\Omega }w=\lambda w+\frac{1}{2}\left(\frac{\mathrm{d}^2w}{\mathrm{d}x^2}k^2w\right)+f(g^2)w+2g^2f^{}(g^2)Re(w).$$
(28)
We have numerically solved (28) by the shooting method with boundary conditions $`w0`$ ($`x\mathrm{}`$) and $`\mathrm{d}w/\mathrm{d}x=0`$ ($`x=0`$). It is clear from (28) that the complex conjugate of the eigenfunction is also an eigenfunction. Namely, if $`w_\mathrm{e}`$ is an eigenfunction of (28) corresponding to the eigenvalue $`\mathrm{\Omega }_\mathrm{e}`$, then $`\overline{w}_\mathrm{e}`$ is also an eigenfunction of (28) corresponding to the eigenvalue $`\overline{\mathrm{\Omega }}_\mathrm{e}`$, where the bar denoting the complex conjugate.

For $`x+\mathrm{}`$, $`g`$ approaches a constant value $`g_0`$, which can be expressed as

$$g_0=\sqrt{1\left(\frac{2}{\lambda +2}\right)^2},$$
(29)
where the inhomogeneity due to the potential vanishes. In this region, $`x`$ also becomes an ignorable direction, and we can assume $`we^{\kappa x}`$ for the point spectra ($`\kappa >0`$). Plugging it into our eigenvalue problem, we obtain the following “dispersion relation” which is applicable for large $`x`$,

$$\mathrm{\Omega }_\mathrm{i}^2=\left[\mathrm{\Omega }_\mathrm{r}+\frac{1}{2}(\kappa ^2k^2)\right]\left[\mathrm{\Omega }_\mathrm{r}+\frac{1}{2}(\kappa ^2k^2)+2g_0^2f^{}(g_0^2)\right],$$
(30)
where the subscripts r and i denote the real and imaginary part, respectively, and we have used the relation $`f(g_0^2)=\lambda `$. By solving it for $`\kappa `$, we obtain

$$\kappa =\sqrt{k^22\mathrm{\Omega }_\mathrm{r}2g_0^2f^{}(g_0^2)\pm 2\sqrt{g_0^4[f^{}(g_0^2)]^2\mathrm{\Omega }_\mathrm{i}^2}},$$
(31)
where $`f^{}(g^2)`$ is explicitly shown as

$$f^{}(u)=(1u)^{3/2}.$$
(32)
The asymptotic argument of $`w`$ is also determined by

$$\mathrm{arg}(w)=\mathrm{arctan}\left(\frac{2\mathrm{\Omega }_\mathrm{i}}{k^2+2\mathrm{\Omega }_\mathrm{r}\kappa ^2}\right).$$
(33)
With these informations, we have solved (28) numerically. First, we fix an arbitrary value on the real part of $`w`$ at $`x=10`$. Given an expectant eigenvalue, the imaginary part of $`w`$ is calculated from (33). Then, the derivative $`\mathrm{d}w/\mathrm{d}x`$ is estimated by using (31), and we can integrate (28) from $`x=10`$ to $`x=0`$ by the 4-th order Runge-Kutta formula. The data of $`g(x)`$ is drawn from the numerical solution of (11). At $`x=0`$, the boundary condition $`\mathrm{d}w/\mathrm{d}x=0`$ is checked. If this boundary condition is not satisfied, we guess the next eigenvalue by the Newton method and carry out the shooting again.

The dispersion relation is shown in Fig. 1. We can see a good qualitative agreement with the analytic evaluation with the fact that the growth rate begins from zero as $`k0`$, experiences a maximum value with respect to $`k`$, and the mode is finally stabilized for a sufficiently large $`k`$. The eigenfunctions corresponding to $`k=0.1`$, $`0.5`$, and $`0.7`$ for the parameter $`\lambda =0.5`$ are illustrated in Fig. 2. From Fig. 1, the mode with $`k=0.5`$ gives the maximum growth rate $`Im(\mathrm{\Omega })0.137`$. As is seen from Fig. 2, the imaginary part of the eigenfunction becomes wider and wider as $`k0`$, which agrees qualitatively with the analytical estimate. Since we are performing the numerical shooting in the finite domain, it becomes difficult to estimate the correct eigenvalue in this regime.

### III.2 Nonlinear evolution

We will now present the numerical simulation of the dynamics of angular momentum creation from one-dimensional dark stripe solitons by solving the full nonlinear system (5)-(6). The calculation was carried out on a $`200\times 200`$ spatial mesh placed in a calculation box of size $`L=28`$. The boundary conditions $`_xA=0`$ ($`_yA=0`$) are imposed at the edges $`x=\pm L`$ ($`y=\pm L`$). Note that these boundary conditions numerically respect the constants of motion in the domain of integration. As an input we choose the one-dimensional $`x`$ dependent stationary solution with $`\lambda =0.5`$. It is to be noted that the photon number and the Hamiltonian are conserved to order $`10^5`$ during the calculation. Since the initial value of the integrated angular momentum is zero, it is found to always remain zero (within an error of order $`10^{15}`$).

Strictly speaking, the eigenfunction in the previous section (see Fig. 2) contains cores of angular momenta, i.e. crossings of the real and imaginary zeros of the field due to the periodic form of the perturbation in the $`y`$ direction. Thus, the growth of small amplitude perturbation itself does not exactly mean the ‘creation’ of angular momenta. In order to check the real creation of angular momenta from an exact zero everywhere, we have first carried out the calculation with an initial condition of the form

$$A(x,y,0)=\alpha \left[g(x)+\stackrel{~}{A}_1\frac{\mathrm{d}g}{\mathrm{d}x}\mathrm{exp}(\beta y^2)\right]+\mathrm{i}\sqrt{1\alpha ^2}g(x),$$
(34)
with $`\alpha =0.9`$, $`\beta =0.0875`$, and $`\stackrel{~}{A}_1=0.1`$. This initial condition is so arranged that there is no crossing of the real and imaginary zeros. It is done by imposing a real-valued perturbation on a complex-valued stationary solution. The zero line of the imaginary part exactly coincides with the $`y`$-axis, while that of the real part deviates from the $`y`$-axis by $`\stackrel{~}{A}_1`$. As the system evolves, we observe the creation of crossings and the appearance of vortex solitons. As the initial conditions are changed (always starting form no crossings), the main qualitative result — the emergence of angular momentum carrying vortex soliton chains — is found to be fairly robust.

Next we investigate the evolution of the system under transverse perturbations of the form

$$A(x,y,0)=g(x)+\stackrel{~}{A}\frac{\mathrm{d}g}{\mathrm{d}x}\mathrm{cos}\left(\frac{n\pi }{L}y\right),$$
(35)
with $`\stackrel{~}{A}=0.1`$ and $`n=6`$. The zero lines of the real and imaginary parts at $`t=2`$ are shown in Fig. 3. As we noted in Sec. II, the crossing points of two zero lines correspond to the vortex centers. We initially have twelve crossings suggesting twelve vortices. These twelve vortices, however, do not have well-formed solitonic structures since they are too close, and overlapping.

As the system evolves, the vortices move and we see annihilations of pairs of opposite polarity. The first annihilation event was observed in the interval $`28\text{ }\stackrel{<}{}\text{ }t\text{ }\stackrel{<}{}\text{ }34`$; two of the vortices destroy one another near the center, and the other two seem to disappear near the edge; the annihilation of two pairs removes four from the original twelve at $`t=2`$.

The second annihilation event occurs during $`40\text{ }\stackrel{<}{}\text{ }t\text{ }\stackrel{<}{}\text{ }48`$. As is depicted in Fig. 4, the zero lines of the real and imaginary parts tend to separate around the center $`y0`$. Here two vortices approach the origin around $`40\text{ }\stackrel{<}{}\text{ }t\text{ }\stackrel{<}{}\text{ }44`$, and then they are annihilated. After the annihilation, two inner vortices approach the origin. Finally the remaining six vortices tend to spread and align with equal inter-vortex spacing (see $`t=72`$ in Fig. 4).

The amplitude $`|A|`$ corresponding to this sequence is illustrated in Fig. 5. At $`t=48`$, the vortices do not quite look like vortex solitons, even though the number of vortices has been reduced to six. Since the distance between adjacent vortices is still small, they overlap and are not quite independent. The central hump around $`xy0`$ is the remnant of the second annihilation event. With time, the peaks tend to expand to the central region and are sufficiently apart to look and behave like vortex solitons. By $`t72`$, we observe the formation of six solitonic structures as predicted in Ref.TBM-01  (see $`t=72`$ in Fig. 5). Notice that the vortices disappear only when annihilated by another with opposite polarity.

After $`t\text{ }\stackrel{>}{}\text{ }80`$, a propagating wave is clearly seen in the region $`y<10`$ and $`y>10`$ where the field was originally flat (see Fig. 6). This propagating field is the trace of the Cherenkov radiation which comes from the non-integrability of the system Iv-97 . The radiation will propagate away from vortices. In the finite calculation domain, the radiation will be reflected back at the boundary.

When we extend the time evolution further, we observe the third annihilation event around $`t160`$. However, this annihilation may be an artifact of the finite size of the domain. When we carry out the simulation in a domain larger than ($`L=28`$) but with the same initial conditions, the third annihilation event takes place at a later time while the times of the first and the second annihilation events remain unchanged. The third annihilation event is likely to be driven by the reflection of the Cherenkov radiation at the boundary. Equivalently, it may be concluded that six vortex solitons are stable in our domain and may stay forever. It is noted that the distance among vortices at the final stage approximately coincides with the inverse of the wave number with the maximum growth rate $`k_\mathrm{m}`$.

## IV Summary

We have demonstrated the dynamics of angular momenta creation in a highly relativistic electron-positron plasmas subject to the passage of a strong pulse of electromagnetic fields. The system is governed by a generalized nonlinear Schrödinger equation with a defocusing inverse square root type nonlinearity. It turns out that the one-dimensional dark soliton stripe solutions of this equation, just like those of the standard nonlinear Schrödinger equation, are unstable to transverse perturbations. By carrying out a linear analysis, we have found that there exists an instability window of transverse wave numbers for the system. By a numerical simulation of the fully nonlinear equation, we have shown that the transverse instability will yield, after a few annihilation events, a well-separated chain of vortex solitons with alternating, singly-charged polarity or topological charge ($`m=\pm 1`$). These singly-charged vortex solitons are topologically stable and do not disappear unless they collide with their compliments and annihilate. The number of the created vortex solitons seem to be determined by the inverse of the wave number $`k_\mathrm{m}`$ with the maximum linear growth rate. For the box size $`L=28`$, the six vortex soliton state is found to be robust.

We have suggested a simple and plausible mechanism of angular momentum generation in the MeV epoch of the Universe. Electromagnetism, operating through the versatile substrate of the electron-positron plasma, seems to readily generate highly interesting, long-lived objects which are capable of carrying large amounts of mass, energy and angular momentum. Since an initial localization of mass, energy and angular momentum is precisely the seed that gravity needs for eventual structure-formation, electromagnetism may have provided a key element in the construction of the large-scale map of the observable Universe.

###### Acknowledgements.

This work is supported in part by a Grant-in-Aid from the Japanese Ministry of Education, Culture, Sports, Science and Technology, No. 12780353. The work of V.I.B. was partially supported by the ISTC grant G-663, and S.M.M.’s work was supported by the U.S. Department of Energy Contract No. DE-FG03-96ER-54346.

Figure Captions

1. Dispersion relation for $`\lambda =0.2`$, $`0.5`$, and $`1`$.
2. Eigenfunctions for $`\lambda =0.5`$ corresponding to $`k=0.1`$, $`0.5`$, and $`0.7`$.
3. Zero lines of real and imaginary parts of the field at time $`t=2`$.
4. Crossing of zeros at time $`t=40`$, 44, 48 and 72.
5. Amplitude $`|A|`$ at time $`t=48`$ and 72.
6. Radiation appears after $`t>80`$.