text:
# Unweighted event generation in hadronic 𝑊𝑍 production at order(𝛼_S)

## I INTRODUCTION

Next-to-leading order in $`\alpha _\mathrm{S}`$ (NLO($`\alpha _\mathrm{S}`$)) corrections in di-boson production are large at the CERN LHC and sizeable at the Fermilab Tevatron (see for example ). Run I physics analyses (i.e., ) at the Tevatron have employed constant $`k`$-factors to approximate NLO($`\alpha _\mathrm{S}`$) corrections. With the increased energy and luminosity of Run II and the new energy regime that will be probed with the CERN LHC, simulations incorporating higher order corrections will be of increased importance.

NLO($`\alpha _\mathrm{S}`$) differential cross-section predictions for di-boson production at hadron colliders are normally accomplished with Monte Carlo integration programs. Divergences inherent in the Feynman graphs contributing at $`𝒪`$($`\alpha _\mathrm{S}`$) are handled with either the phase space slicing method  or the subtraction method  with the effect that both positive weighted events and negative weighted (or probability) events are generated. The contribution from negative weighted events cancels in an integration over a suitably large volume of hadronic final states, and a physical (positive probability) prediction is obtained. Integration programs are successful in producing distribution and cross-section predictions for inclusive observables subject to any experimental cuts providing statistics are large enough to effect the cancellations, and a sufficient volume of jets has been integrated within each histogram bin. However, owing to the negative probability events, unweighted event generation — the generation of events with the distribution predicted by theory — with these NLO($`\alpha _\mathrm{S}`$) matrix elements has not been achieved to date.

Unweighted event generation is often preferred over integrated distributions because it can genuinely simulate experimental data and because it is difficult to perform detector simulation on the large event statistics required for integrated distributions.

Unweighted event generation can be reduced to a problem of rendering all event weights positive definite. If a sample of positive definite weighted events can be obtained, the Monte Carlo hit-and-miss method is commonly used for unweighting. The maximum event weight $`d\sigma _{\mathrm{m}ax}`$ for the process is estimated, for instance, by sampling the cross-section a number of times. The ratio of event weight over the maximum event weight $`\frac{d\sigma }{d\sigma _{\mathrm{m}ax}}`$ is compared to a random number $`g`$ generated uniformly in the interval (0,1). Events for which the ratio exceeds the random number ($`\frac{d\sigma }{d\sigma _{\mathrm{m}ax}}>g`$) are accepted, the others are rejected. This ensures the accepted events are distributed according to the description provided by the matrix element, and so all have unit weight.

In this paper we describe an algorithm which extends the usefulness of the present NLO($`\alpha _\mathrm{S}`$) hadronic di-boson production programs by providing a means for unweighted event generation thereby presenting the NLO($`\alpha _\mathrm{S}`$) predictions in a format better suited for experimental simulations. It deals specifically with the region of low transverse momentum of the $`WZ`$ system, where higher order corrections are largest. As for the Drell-Yan process, resummation (first proposed by ) is expected to give more reliable predictions in this region — however a resummed treatment of the $`WZ(j)`$ process does not exist and a comparison with fixed order predictions will always be of interest. Another approach for this problematic region is the matching of the zero-order matrix elements using a parton shower model for simulations of the initial-state shower to the first order tree level (i.e., $`𝒪`$($`\alpha _\mathrm{S}`$)) matrix elements. This has been accomplished for the single vector boson case .

The algorithm makes use of the matrix elements from the Baur, Han, and Ohnemus (BHO)  integration package which employs the phase space slicing method as regularization scheme. This scheme makes use of approximations for the analytic integration of a volume of soft and collinear real emissions. The approximations break down before the volume can be increased sufficiently to render all event weights positive.

The algorithm presented here is an extension of the phase space slicing method. A cylinder of jets defined by a transverse momentum cutoff $`P_{\mathrm{c}utoff}^T`$ is used to partition event generation into a subprocess where the real emission has sufficient transverse momentum to be potentially observable (called 1 jet) and another subprocess (called 0 jet) where the real emission (if present) has sufficiently small transverse momentum such that it is taken as unobservable. There are no divergences in the 1 jet cross-section so these events all have positive weights and unweighted event generation may be accomplished by applying Monte Carlo hit-and-miss. Difficulties arise in the treatment of the 0 jet events, which is the topic of this paper.

The missing transverse momentum in 0 jet events may arise from the invisible neutrino and from jets which escape detection. By fixing the total missing transverse momentum rather than the individual contributions from the invisible neutrino and unobserved jets, an integration volume of real emissions phase space can be generated and so the cancellations between the real emissions and one-loop graphs can be accomplished by a Monte Carlo integration over this jet volume. The charged lepton and missing transverse momentum vectors (i.e., the observable quantities) are kept fixed in the laboratory frame throughout the integration. Since the jet volume may be large, the fixing of the observable quantities is essential. An integration which instead keeps other quantities such as the subprocess energy and boost fixed (as would be the case if the phase space slicing method were applicable to such large jet volumes) would also necessarily be an integration over observationally different charged lepton vectors.

A final state neutrino is not the only means of generating the missing transverse momentum over which a jet volume can be integrated. Anything which contributes to the missing transverse momentum may be used including the primordial $`k^T`$ distribution of the proton, multiple interactions, minimum bias events, and even detector resolution. For $`WZ`$ production with leptonic decays, the invisible neutrino is a convenient choice. In generalizing this algorithm to other processes different contributions to the missing transverse momentum could be used.

Section II of this paper briefly enumerates the presently available integration packages for hadronic $`WZ`$ production and outlines the regularization schemes commonly used for the matrix elements before describing the algorithm in some detail. In Section III cross-sections and several differential cross-section distributions from integration packages are compared to events generated with this algorithm. Our conclusions are presented in the final section.

## II ALGORITHM FOR UNWEIGHTED EVENT GENERATION

We choose the specific process $`ppWZ(j)`$ as a case study to apply the algorithm. LHC energy is chosen because the challenges in terms of negative probability events are more severe at 14 TeV than for the Tevatron case.

### A Order $`\alpha _\mathrm{S}`$ cross-section

Hadronic $`WZ`$ production at NLO($`\alpha _\mathrm{S}`$) has been calculated in Ref. . All calculations with the exception of Ref.  use the narrow width approximation, wherein the gauge bosons are taken on their mass shell. The correlations in the decays of the vector bosons to massless fermions are fully included in Ref. , while they are included in all but the one-loop graphs in Ref. . For the numerical implementation of the NLO($`\alpha _\mathrm{S}`$) matrix elements Ref.  uses the subtraction method, while Ref.  employs the phase space slicing method. The calculations of Ref.  were found to be in agreement in Ref. . An updated version of the Ref.  program incorporating an omitted region of phase space has been compared with Ref.  and found to be in good agreement \[1, Sec. 5.5\].

The NLO($`\alpha _\mathrm{S}`$) cross-section receives contributions from the square of the Born level graphs, the interference of the Born graphs with the one-loop graphs, and the square of the real emission graphs,

$$_{\mathrm{N}LO}^2=_{\mathrm{B}orn}^2+_{\mathrm{B}orn}_{\mathrm{O}neLoop}+_{\mathrm{R}ealEmission}^2.$$
(1)
Soft, collinear, and ultraviolet divergences appear when any of the $`𝒪(\alpha _\mathrm{S})`$ graphs are treated alone, and so a regularization scheme is necessary to ensure the inherent cancellations. The matrix elements used in this study employ the phase space slicing method. A volume of real emissions phase space (herein referred to as n+1 body phase space) is partitioned by means of soft and collinear cutoffs The phase space slicing method defines the soft jet region by
$$E_{\mathrm{j}et}^{\mathrm{C}MS}<\delta _\mathrm{s}\frac{\sqrt{\widehat{s}}}{2}$$
where $`E_{\mathrm{j}et}^{\mathrm{C}MS}`$ is the jet energy in the $`WZj`$ subprocess center of mass frame, $`\delta _\mathrm{s}`$ is the soft cutoff, and $`\sqrt{\widehat{s}}`$ is the subprocess energy. The collinear region is defined by
$$|E_iE_j\pm \stackrel{}{p}_i\stackrel{}{p}_j|<\delta _\mathrm{c}\widehat{s}$$
where $`E_i`$ and $`\stackrel{}{p}_i`$ are the energy and momentum of massless parton $`i`$ and $`\delta _\mathrm{c}`$ is the collinear cutoff. For more details see Ref. .  and incorporated analytically with the soft gluon and leading pole approximations into the phase space of the diagrams lacking a real emission (herein referred to as n body).

In the region where the soft and collinear divergences overlap, the perturbative QCD calculation may be visualized in terms of an expansion in powers of

$$\frac{\alpha _\mathrm{S}}{2\pi }\left[\mathrm{ln}\frac{\widehat{s}}{P_{\mathrm{j}et}^{T}{}_{}{}^{2}}\right]^2$$
(2)
where $`\alpha _\mathrm{S}`$ is the QCD coupling, $`\sqrt{\widehat{s}}`$ is the subprocess energy and characterizes the scale of the interaction, and $`P_{\mathrm{j}et}^T`$ is the transverse momentum of the order $`\alpha _\mathrm{S}`$ real emission. When the jet transverse momentum $`P_{\mathrm{j}et}^T`$ is small relative to $`\sqrt{\widehat{s}}`$, the function given in Eq. 2 becomes large and multi-gluon emission is important — this is the problem region for fixed order calculations. We wish to integrate out a large enough $`P_{\mathrm{j}et}^T`$ volume such that the function given in Eq. 2 remains smaller than unity. For $`\sqrt{\widehat{s}}=1`$ TeV this gives $`P_{\mathrm{j}et}^T`$ of order 10 GeV, which in turn corresponds nicely to the detector capabilities at the LHC, wherein typically jets above 15 (30) GeV will be reconstructed at low (high) luminosity. Smaller jet integration volumes would be necessary at lower machine energies.

The situation is illustrated schematically in Fig. 1. A hypothetical jet transverse momentum distribution is shown for the n body and n+1 body case. The n+1 body contribution diverges as the jets become soft or collinear, while the n body contribution is a negative divergence at the origin. The divergences cancel in an integration over a suitably large volume of jet transverse momenta (hatched region). The phase space slicing method effects such an integration, but the jet volume is not sufficiently large to keep the cross-section in the inclusive (hatched) region positive for most events.

### B Generating unweighted events

We begin by dividing events into two classes: those containing a jet with sufficient transverse momentum ($`P_{\mathrm{j}et}^T>P_{\mathrm{c}utoff}^T`$) to be classified as potentially observable (referred to here as 1 jet); and those without a jet or containing a jet with transverse momentum below the cutoff $`P_{\mathrm{j}et}^T<P_{\mathrm{c}utoff}^T`$ (referred to as 0 jet). This differs from the definition of n body and n+1 body events where the soft and collinear cutoffs define the partition. In the 1 jet region, the matrix elements are order $`\alpha _\mathrm{S}`$, but the calculation is leading order. There are no divergences nor are there negative weighted events, so event generation can be accomplished with the usual Monte Carlo techniques.

The algorithm concerns events in the 0 jet region only and we focus exclusively on events of this type for the remainder of this paper.

The phase space slicing method is extended by integrating out the cylinder of jets defined by the $`P_{\mathrm{c}utoff}^T`$ using a second stage Monte Carlo integration. As the $`P_{\mathrm{c}utoff}^T`$ is chosen to be the order of 10 GeV, there is a danger that the integration may affect the kinematics of the event to the point which is observable by the detector. With this in mind, the observable vectors (charged leptons for the case of $`WZ`$ production with leptonic decays) are kept fixed throughout the integration, and the center of mass frame vectors are allowed to fluctuate.

We start by sampling the observable quantities of the event i.e., fixing the vectors

$$P_{\mathrm{l}_W^\pm },P_{\mathrm{l}_Z^+},P_{\mathrm{l}_Z^{}},\stackrel{}{P}_{\mathrm{m}iss}^T$$
(3)
where $`P_{\mathrm{l}_W^\pm }`$ is the momentum of the charged lepton from the $`W^\pm `$ decay, $`P_{\mathrm{l}_Z^+},P_{\mathrm{l}_Z^{}}`$ are the momenta of the charged leptons from the $`Z`$ decay, and $`\stackrel{}{P}_{\mathrm{m}iss}^T`$ is the missing transverse momentum. With the $`Z`$ boson taken on shell and requiring no net momentum transverse to the beam, this amounts to 8 degrees of freedom.

For the case of the n body matrix element these vectors also fix the neutrino four-vector with $`\stackrel{}{P}_\nu ^T=\stackrel{}{P}_{\mathrm{m}iss}^T`$ and the $`W`$-mass constraint gives zero or two solutions for the neutrino longitudinal momentum. Within the BHO matrix elements a further 4 degrees of freedom specify the NLO($`\alpha _\mathrm{S}`$) corrections. Thus the average n body weight $`d\sigma ^{\mathrm{n}\mathrm{body}}`$ for the observable event is arrived at by integrating over these degrees of freedom and averaging over the two neutrino solutions. This weight is usually negative. When no neutrino solutions exist the weight is zero. We refer to a point in the full n body phase space (defined by the vectors of Eq. 3, the 4 extra degrees of freedom, and the neutrino solution choice) as a “subevent”, such that an “event” as defined by this algorithm is an integration over many subevents.§§§ A subevent from this algorithm corresponds precisely to an event from a Monte Carlo integration package.

We now incorporate the cylinder of jets defined by the $`P_{\mathrm{c}utoff}^T`$ by performing a similar integration over n+1 body subevents. For each subevent a jet with transverse momentum up to the cutoff is sampled (3 degrees of freedom). The jet transverse momentum fixes the invisible neutrino transverse momentum for the event with the constraint that the observed missing transverse momentum is kept fixed,

$$\stackrel{}{P}_\nu ^T=\stackrel{}{P}_{\mathrm{m}iss}^T\stackrel{}{P}_{\mathrm{j}et}^T.$$
(4)
The longitudinal neutrino momentum is specified up to a two-fold ambiguity by the $`W`$-mass constraint. The charged lepton vectors are also kept fixed in the laboratory frame, such that all the quantities in Eq. 3 are unchanged throughout the integration. The average subevent weight is $`d\sigma ^{\mathrm{n}+1\mathrm{body}}`$, which is a positive number.

The phase space for a Monte Carlo integration is normally sampled in terms of variables such as the invariant masses, rapidities, and decay angles which give a triangular transformation matrix such that the Jacobian is the product of the diagonal entries (as is the case for the di-boson integration packages). The fixed charged lepton vectors in the laboratory frame constraint introduces inherent correlations, and the transformation matrix becomes more complicated: for the $`WZ`$ 1 jet case it is a $`11\times 11`$ matrix with no zero elements. Fortunately it is simple to calculate this matrix and its Jacobian numerically on an event by event basis, which is the strategy we have used. This Jacobian is included in the subevent weights $`d\sigma ^{\mathrm{n}\mathrm{body}}`$ and $`d\sigma ^{\mathrm{n}+1\mathrm{body}}`$.

The total event weight for the configuration defined by Eq. 3 is the sum of the n body and n+1 body contributions, with a factor which accounts for the observable vectors phase space sampling

$$d\sigma ^{0\mathrm{jet}}=W_{\mathrm{P}haseSpace}\times \left(d\sigma ^{\mathrm{n}\mathrm{body}}+d\sigma ^{\mathrm{n}+1\mathrm{body}}\right).$$
(5)
It is unimportant that the number of subevents for the second stage integration be sufficiently large to give an accurate result for each observable event. A poor accuracy simply means two identically configured events may have differing event weights, however by grace of the Monte Carlo method, the average converges to give the correct cross-section. The subevent statistics are required to be sufficiently large to ensure the event weights of Eq. 5 very rarely fall negative on account of statistical fluctuations. We refer to events where this occurs as “remnant negative weight events”. Events of this type are discarded, and thus contribute to a (biased) systematic Monte Carlo error. The number of subevents is a means of directly controlling this error. It is easy to evaluate this error.

Finally, we note the application of this algorithm to the re-weighting of leading order events (i.e., re-scaling of tree level event weights to account for NLO($`\alpha _\mathrm{S}`$) corrections) needs to be approached with caution, as the phase space volume of observable vectors in Eq. 3 is larger for the n+1 body case than for the n body case since only limited configurations of these vectors provide n body neutrino longitudinal momentum solutions which satisfy the $`W`$-mass constraint. This difference needs to be accounted for in a re-weighting scheme.

## III RESULTS AND COMPARISON WITH INTEGRATION PACKAGES

The results presented in this section use identical parameters as the comparison of integration packages in Ref. \[1, Sec. 5.5\] with the exception of the jet definition for the jet veto. In that study the jet transverse momentum cutoff $`P_{\mathrm{c}utoff}^T`$ was motivated by what the detectors are capable of measuring. Here we motivate the $`P_{\mathrm{c}utoff}^T`$ by what the detector is incapable of reconstructing: a $`P_{\mathrm{c}utoff}^T`$ of 15 GeV is defined to partition the 1 jet and 0 jet regions (there is no jet rapidity requirement). The CTEQ4M  structure functions are used. Input parameters are taken as $`\alpha _{EM}(M_Z)=\frac{1}{128}`$, $`\mathrm{sin}^2\theta _W=0.23`$, $`\alpha _s(M_Z)=0.1116`$, $`M_W=80.396`$ GeV, $`M_Z=91.187`$ GeV, factorisation scale $`Q^2=M_W^2`$, and Cabibbo angle $`\mathrm{cos}\theta _{\mathrm{C}abibbo}=0.975`$ with no 3rd generation mixing. Branching ratios are taken as Br($`Zl^+l^{}`$) = 3.36%, Br($`W^\pm l^\pm \nu `$)=10.8%. The $`b`$ quark contribution to parton distributions has been taken as zero. Kinematic cuts motivated by anomalous triple gauge boson coupling analyses are chosen (see for example \[1, Sec. 5\]). The transverse momentum of all charged leptons must exceed 25 GeV and the rapidity of all charged leptons must be less than 3 in magnitude. Missing transverse momentum must be greater than 25 GeV. The BHO package is employed for the evaluation of matrix elements.

We first evaluate the algorithm’s ability to integrate out the negative probability events. The second stage integration is accomplished with relatively low statistics to keep the computer processing time reasonable, only 100 subevents of each type (n body and n+1 body) are used.

For the default choices of the soft and collinear cutoffs in the BHO matrix elements ($`\delta _\mathrm{s}=0.01,\delta _\mathrm{c}=0.001`$) the cancellations between the n body and n+1 body partitions are very large, and the second stage integration would need unreasonably large statistics for good behavior.The same problematic behavior is observed when matrix elements employing the subtraction method are used. We found no immediate solution to this problem for the subtraction method, and so we have been unable to use this algorithm with subtraction method matrix elements thus far. In order to achieve a portion of the jet volume integration analytically we increase the cutoffs ($`\delta _\mathrm{s}=0.05,\delta _\mathrm{c}=0.002`$) which is well within the range of applicability of the soft gluon and leading pole approximations. We emphasize that experiments must take the integration volume into account by ensuring that their analysis is not sensitive to jets within the volume (as defined by the jet momentum cutoff of this algorithm and the soft and collinear cutoffs of the phase space slicing method). This consideration is true in general for the phase space slicing method.

With these choices of the phase space slicing cutoffs “remnant negative weight events” account for approximately 0.1% of the total sample of weighted events and -0.02% of the cross-section. Most of these events are the result of the low second stage integration statistics such that re-evaluating these events with a larger number of subevents renders the event weight positive.

There is another class of remnant negative events which are not due to statistics. For these events the charged lepton in the $`W`$-decay is emitted along the line of flight of the $`W`$ in its rest frame. This configuration is disfavored by the spin correlations, wherein the $`W^+`$-decay ($`W^{}`$-decay) with a left-handed (right-handed) particle emitted along the $`W`$ line of flight is preferred. The BHO matrix elements use a spin averaged treatment for the one-loop graphs and so these correlations are not at all included for the virtual corrections (though spin correlations are included everywhere else). This means events in this configuration receive an unphysically large negative contribution from the one-loop graphs while receiving the correct (spin correlation suppressed) contribution from the other graphs. Thus this class of negative event weights arise as a result of the lepton correlation approximations used in the BHO matrix elements. The algorithm is sensitive to these approximations and would benefit from improved theoretical calculations (as are included in ). The negative weights associated to events of this type are extremely small.

Finally there are a few events with negative weights at extraordinary subprocess energy ($`\sqrt{\widehat{s}}>5`$ TeV) where the total cross-section is orders of magnitude below one event per LHC year. These events may arise from a combination of the lepton correlation effect and an insufficiently large jet volume for events at this energy scale, as suggested by the function given in Eq. 2.

### A Comparison with the BHO integration package

In Table I we enumerate the leading order and NLO($`\alpha _\mathrm{S}`$) predictions for 0 jet $`WZ`$ production from the BHO integration package together with the result from this algorithm which employs the BHO matrix elements. The NLO($`\alpha _\mathrm{S}`$) prediction differs considerably from the leading order one and this algorithm is able to reproduce the NLO($`\alpha _\mathrm{S}`$) BHO predictions well. The accuracy of the prediction from this algorithm suffers slightly with respect to that from the BHO integration package alone in that the correlations between integration variables make it more difficult for the adaptive integration to optimize the phase space sampling.

It takes considerably more computer time to evaluate the NLO($`\alpha _\mathrm{S}`$) matrix elements than for the leading-order case. This algorithm requires a further factor 200 over the NLO($`\alpha _\mathrm{S}`$) integration package in evaluating one weighted event because of the 2$`\times `$100 subevents in the second stage integration.

Monte Carlo hit-and-miss is used to extract a sample of unweighted events from the (positive definite) events generated using this algorithm. No attempt has been made to optimize for efficiency, and the time required for the generation of one unweighted event is quite large 200 seconds are required for one unweighted event using a 600 MHz personal computer running Linux 2.2. Pythia  requires only about 0.02 seconds for the generation of one leading order $`WZ`$ event with an efficiency of 1-2%.  with an efficiency of about 0.1%. Nevertheless it is possible to generate an event sample of a few low luminosity LHC years (i.e., several thousand events) in several hours using computer farms of the size currently in use at most universities or experiments.

Weighted 0 jet event distributions are compared in Fig. 2 & 3 for the BHO integration package and this algorithm. Figure 2(top) shows the transverse momentum distribution of the $`Z`$, which is the distribution normally used to probe anomalous $`WWZ`$ couplings at hadron colliders . The rapidity separation between the $`Z`$-boson and the charged lepton from the $`W`$-decay (Fig. 2, bottom) is sensitive to the approximate radiation zero in hadronic $`WZ`$ production . Finally the missing transverse momentum distribution is shown in Fig. 3. The jet integration volume arises from this distribution, and so it is the place where discrepancies between this algorithm and the integration package might arise. The agreement between the BHO integration package and this algorithm in all of these distributions is good.

In order to evaluate the effect of remnant negative weight events, the absolute value of the weights of negative events is plotted in Fig. 2 & 3 (shaded regions). The remnant negative weight event distribution is several orders of magnitude below the positive event distribution in all cases and is largest in the region where the differential cross-section is also largest, and so the effect of these events is negligible everywhere. Histograms of this sort provide the means of evaluating the effect of these remnant negative events. For faster event generation the second stage integration statistics may be decreased, resulting in an increase in the size of this shaded region. It is impossible to completely eliminate the shaded region of remnant negative events by increasing the second stage integration statistics indefinitely, due to the lepton correlation effects discussed in the first part of this section. If the negative weighted events from the BHO integration package (without this algorithm) were also superimposed, they would be distributed uniformly about half an order of magnitude below the BHO prediction.

In Fig. 4 unweighted events from this algorithm are superimposed on the BHO integration package prediction for the same distributions as Fig. 2. The error bars assume Poisson statistics and the agreement is good. This figure represents the first unweighted event generation at NLO($`\alpha _\mathrm{S}`$) for the partonic process $`ppWZ(j)`$.

## IV CONCLUSIONS

We have described an algorithm for unweighted event generation at NLO($`\alpha _\mathrm{S}`$) for the partonic process $`ppWZ(j)`$ with leptonic decays and evaluated its effectiveness at LHC energy. Event generation consists of a two stage Monte Carlo integration, and so requires considerable computer time. The possibility of optimizing the generation and the code itself has not been addressed, and certainly the performance may be considerably improved. The partition of events into 0 jet and 1 jet regions is motivated by experimental capabilities and defined by a jet transverse momentum cutoff. This partition, as well as the n body and n+1 body partition inherent in the phase space slicing method, must be taken into account in experimental analyses.

Generation of high transverse momentum real emissions is not addressed by this algorithm, but this region of events is free of divergences and so unweighted event generation is easily achieved with the usual Monte Carlo methods.

A small number of remnant negative weighted events are discarded by the algorithm. The effect of these events is negligible and contributes to a biased systematic Monte Carlo error which may be monitored and evaluated.

The agreement between events generated with this algorithm and the BHO NLO($`\alpha _\mathrm{S}`$) integration package is good.

The algorithm is generalizable to a broader class of processes. Hadronic $`W\gamma (j)`$ production can be treated in an identical way, while the $`W^+W^{}(j)`$ process with leptonic decays can be handled by treating the neutrino (anti-neutrino) as invisible and the anti-neutrino (neutrino) as an observable lepton for the purposes of this algorithm. Other means of generating the jet integration volume besides the invisible neutrino may be used for other processes.

## ACKNOWLEDGMENTS

The authors would like to thank the ATLAS collaboration which motivated this project. We are grateful to Ulrich Baur for providing the matrix elements as well as informative correspondence and discussions. We thank Ian Hinchliffe and J$`ø`$rgen Beck Hansen for useful comments. This work has been supported by the Natural Sciences and Engineering Research Council of Canada.

## FIGURES

## TABLES


format: text

text:
# Angle-resolved photoemission study of insulating and metallic Cu-O chains in PrBa2Cu3O7 and PrBa2Cu4O8

\[

## Abstract

We compare the angle-resolved photoemission spectra of the hole-doped Cu-O chains in PrBa2Cu3O7 (Pr123) and in PrBa2Cu4O8 (Pr124). While, in Pr123, a dispersive feature from the chain takes a band maximum at $`k_b`$ (momentum along the chain) $``$ $`\pi /4`$ and loses its spectral weight around the Fermi level, it reaches the Fermi level at $`k_b`$ $``$ $`\pi /4`$ in Pr124. Although the chains in Pr123 and Pr124 are approximately 1/4-filled, they show contrasting behaviors: While the chains in Pr123 have an instability to charge ordering, those in Pr124 avoid it and show an interesting spectral feature of a metallic coupled-chain system.

\]

One-dimensional (1D) or quasi 1D compounds have attracted great interest because of their unique and interesting physical properties such as Peierls instability, spin-charge separation, and Tomonaga-Luttinger (TL) liquid behavior . In order to realize the TL liquid state, the system should have good one dimensionality. On the other hand, a good 1D system is expected to have instability towards spin and/or charge ordering and tends to be insulating. In every 1D system, a metallic state which is possibly a TL liquid, is competing with the instability to the insulating state. It is therefore very important and interesting to study the electronic structure of quasi 1D systems which are close to the boundary between the quasi 1D metal and insulator. This subject is also related to the electronic structure of the quasi 1D stripe phase in 2D systems and is important in the light of the competition between various quantum mechanical ground states in 1D and 2D systems .

Angle-resolved photoemission spectroscopy (ARPES) is a powerful tool to study spin-charge separation as well as band gap opening due to charge ordering. It has been found that ARPES of the Cu-O double chain in SrCuO2 and the Cu-O single chain in Sr2CuO3 have two dispersive features which are successfully interpreted as spinon and holon dispersions of the undoped Cu-O chains and are manifestation of the spin-charge separation in 1D Mott insulators . In an ARPES study of PrBa2Cu3O7 (Pr123), it has been reported that the hole-doped Cu-O single chain in Pr123 shows possible spinon and holon dispersions and also has band gap opening at Fermi level ($`E_F`$) probably due to charge ordering . Recently, an ARPES study of a quasi 1D metal Li0.9Mo6O17 has shown that spectral weight near $`E_F`$ is considerably suppressed and can be interpreted as a TL liquid with large $`\alpha `$ which is the anomalous exponent of the momentum distribution function as well as of the single-particle density of states . The large $`\alpha `$ suggests strong fluctuations of charge ordering or charge density wave (CDW) in Li0.9Mo6O17 although the analysis of the ARPES data is still controversial .

In Pr123 and PrBa2Cu4O8 (Pr124), the CuO2 planes remain antiferromagntic and insulating and do not have enough carriers to cause superconductivity . On the other hand, the Cu-O single chains in Pr123 and the Cu-O double chains in Pr124 are heavily hole-doped and show semiconducting and metallic behaviors, respectively . Therefore, Pr123 and Pr124 give us a unique opportunity to study the electronic structures of the hole-doped Cu-O chains which are close to the boundary between the quasi 1D metal and insulator. In this Letter, we report on a new set of ARPES data for Pr124 which has metallic Cu-O double chains. The ARPES data of the Cu-O double chain in Pr124 shows a sharp contrast to those of the Cu-O single chain in Pr123. Since the structures of the Cu-O chains in Pr123 and Pr124, which are the same as the Cu-O chains in Sr2CuO3 and SrCuO2 respectively, are simple compared to Li0.9Mo6O17, the interpretation of the ARPES data is rather straightforward. By comparing the ARPES data of Pr123 and Pr124, we discuss the effect of interchain coupling and the instability towards charge ordering in the hole-doped Cu-O chains.

Naturally-untwinned single crystals of Pr124 were grown by a flux method under oxygen pressure of 11 atm. The resistivity along the chain direction is metallic as reported in the literature . The ARPES measurements of Pr124 were performed at beamline 5-4 equipped with a Scienta SES 200 electron analyzer, Stanford Synchrotron Radiation Laboratory (SSRL). The chamber pressure during the measurements was less than 5 x 10-11 Torr. The samples were cooled to 10 K and cleaved in situ. The cleaved surfaces were the $`ab`$-plane, where the $`b`$-axis is in the Cu-O chain direction. The cleanliness of the surfaces was checked by the absence of a hump at $``$ 9.5 eV. The position of $`E_F`$ was calibrated with gold spectra. The experimental uncertainty in the energy calibration was $`\pm `$ 1 meV. The ARPES data of Pr123 were taken at beamlines 5-3 of SSRL. The details of the measurements were described in the previous paper . For the ARPES data of Pr123 and Pr124 shown in this paper, incident photons were linearly polarized and had an energy of 29 eV for Pr123 and 22.4 eV for Pr124. The total energy resolution including the monochromator and the analyzer was approximately 40 meV for Pr123 and 20 meV for Pr124. The angular resolution was $`\pm 1`$ degree for Pr123 and $`\pm 0.28`$ degree for Pr124, which gives the momentum resolution of $`\pm `$ 0.05$`\pi `$ for Pr123 and $`\pm `$ 0.01$`\pi `$ for Pr124 in units of $`1/a`$ or $`1/b`$. ($`a`$ = 3.87 $`\mathrm{\AA }`$ and $`b`$ = 3.93 $`\mathrm{\AA }`$ for Pr123, $`a`$ = 3.88 $`\mathrm{\AA }`$ and $`b`$ = 3.90 $`\mathrm{\AA }`$ for Pr124)

The experimental arrangement is schematically shown in Fig. 1. The polarization vector of the incident photons had a component parallel to the Cu-O chain direction. The CuO4 square planes of the Cu-O chains are perpendicular to the cleaved surface. In Fig. 2, the ARPES spectra of Pr124 along the Cu-O chain direction are shown for $`k_a=\pi `$ . Here, $`k_a`$ and $`k_b`$ are the momentum perpendicular to the chain in units of $`1/a`$ and the momentum along the chain in units of $`1/b`$, respectively. In Fig. 2, one can see that a dispersive feature from the Cu-O chain moves to $`E_F`$ in going from $`k_b/\pi =0.1`$ to $`k_b/\pi =0.2`$. For $`0.2<k_b/\pi <0.25`$, the dispersion becomes flatter and the dispersive feature gradually loses its intensity. This behavior is clearly seen in the right panel of Fig. 2, where the spectral weight integrated from -0.05 eV to $`E_F`$ is plotted as a function of $`k_b/\pi `$. The integrated spectral weight decreases from $`k_b/\pi =0.20`$ to 0.27, indicating that the dispersive feature crosses $`E_F`$ at $`k_b/\pi =0.23\pm 0.03`$  and that the hole concentration of the Cu-O chain in Pr124 is 0.46 $`\pm `$ 0.06. This is consistent with the recent optical study which shows that the hole concentration of the Cu-O chain in Pr124 is $``$ 0.4 .

The ARPES spectra of Pr123 and Pr124 along the Cu-O chain are compared in Fig. 3 for $`k_a=\pi `$. In Pr123, the dispersive feature from the Cu-O chain reaches the band maximum of -0.3 eV at $`k_b/\pi =0.24`$ and stays there at $`k_b/\pi =0.29`$. This feature gradually loses its weight for $`k_b/\pi >0.29`$ without reaching $`E_F`$. On the other hand, in Pr124, the dispersive feature reaches $`E_F`$ at $`k_b/\pi 0.23\pm 0.03`$. Although the Cu-O chains are approximately 1/4-filled both in Pr123 and in Pr124, the spectral feature near $`E_F`$ of Pr124 is remarkably different from that of Pr123. As shown in Fig. 4, the dispersions of the Cu-O chain features become clear in the density plots of the raw ARPES spectra. While, in Pr124, the width of the Cu-O band below $`E_F`$ is $``$ 0.5 eV and the band dispersion is relatively large, the band width is reduced to $``$ 0.3 eV in Pr123 because of the band gap opening at $`E_F`$.

For Pr123, the band gap opening can be attributed to the charge ordering or CDW formation in the 1/4-filled Cu-O single chain . Actually, charge instability in the Cu-O chain has been observed by NMR and NQR measurements of Pr123 . The CDW picture is also consistent with the experimental result that the gap opening at $`E_F`$ reduces the band width of Pr123 compared to Pr124. On the other hand, the Cu-O double chain of Pr124 does not have charge ordering and is metallic. A small deviation from the 1/4-filling would be responsible for this suppression of charge ordering in Pr124. It is also possible that, since each double chain of Pr124 consists of two single chains (see Fig. 1), charge ordering is unstable because of the weak interaction between the two single chains in each double chain .

The ARPES study of Li0.9Mo6O17 shows that the spectral weight near $`E_F`$ is considerably suppressed even in the metallic region and gives the anomalous exponent $`\alpha `$ of 0.9 , indicating that Li0.9Mo6O17 is a TL liquid close to the CDW instability. On the other hand, the ARPES data of Pr124 have substantial spectral weight at $`E_F`$ compared to Li0.9Mo6O17  and, in the low-energy region, Pr124 does not show the TL-liquid behavior. This is consistent with the fact that Pr124 has large Hall coefficient at low temperature and behaves as a 2D system . Probably, the deviation from the TL-liquid behavior is caused by the hopping between the double chains which becomes relevant in the low-energy region.

The hopping term between the double chains is estimated to be $``$ 10 meV from the 1D-2D crossover observed in the transport measurements of Pr124,. In the energy region higher than the hopping term, the 1D character is expected to manifest in the ARPES spectra. In order to show the dispersions more clearly, the second derivatives of the ARPES spectra are displayed in Fig. 5. In Pr123, two dispersive features labeled as $`\alpha ^{}`$ and $`\alpha ^{\prime \prime }`$ are visible as two bright belts which can be attributed to holon and spinon dispersions, respectively, and are the manifestations of the spin-charge separation in 1D systems as predicted theoretically . The holon and spinon dispersions have the width of $``$ 0.5 and $``$ 0.1 eV, respectively, which approximately agree with $`t`$ and $`J`$ in the $`t`$-$`J`$ model for the cuprates . On the other hand, Pr124 does not show separate spinon and holon features expected for the 1D $`t`$-$`J`$ model . One possible explanation is that, in the high-energy region, Pr124 behaves as a TL-liquid with anomalous exponent $`\alpha `$ larger than 0.5 because of proximity to charge ordering as proposed in the recent optical study on Pr124 . In this case, since the TL model predicts that the spinon feature appears as a cusp instead of a power-law divergence, the spinon feature is expected to be broad and may not separately be observed from the holon feature .

As seen in Figs. 3 and 4. the ARPES spectra for $`k_b/\pi `$ $`<`$ 0.20 have substantial spectral weight near $`E_F`$ although the dispersive feature is located well below $`E_F`$ (see Fig. 2). In the second derivative shown in Fig. 5, this spectral weight near $`E_F`$ appears as a horizontal bright belt near $`E_F`$ ranging from $`k_b/\pi `$ = 0.0 to 0.4. In addition, for $`0.20<k_b/\pi <0.25`$, the dispersion becomes flatter and it looks like that a small pseudo-gap tend to open at $`E_F`$. These behaviors in the metallic Cu-O chains of Pr124 are similar to those observed along $`(\pi ,0)(\pi ,\pi )`$ above $`T_c`$ in the underdoped CuO2 plane of the high-$`T_c`$ cuprates , which in the stripe phase represent the dispersion along the 1/4-filled stripe . The similarity between the coupled chains and the high-$`T_c`$ cuprates has been argued by Kopietz, Meden, and Schönhammer in the light of the crossover between a TL liquid and a 2D Fermi liquid . It would be interesting to investigate whether the coupled-chain model can explain the present ARPES data of Pr124 as well as those of the high-$`T_c`$ cuprates.

In conclusion, we have studied the hole-doped Cu-O chains in Pr123 and Pr124 using ARPES. The ARPES data show that the Cu-O single chain in Pr123 and the Cu-O double chain in Pr124 are approximately 1/4-filled. While, in Pr123, the 1D features from the Cu-O single chain lose their weight near $`E_F`$, the dispersive feature from the Cu-O double chain reaches $`E_F`$ in Pr124. These facts indicate that the charge ordering occurs in the Cu-O single chain and is suppressed in the Cu-O double chain. In the low-energy region, the line shape of the ARPES spectra of Pr124 cannot be explained neither by the TL-liquid picture nor by the Fermi-liquid picture and looks rather close to that of the 2D CuO2 plane with stripes. It would be useful to further study the differences and similarities between the metallic Cu-O chains in Pr124 and the metallic stripe phase in the high-$`T_c`$ cuprates.

The authors would like to thank K. Takenaka, I. Terasaki, K. Penc, T. Thoyama, and S. Maekawa for valuable comments and the staff of SSRL for technical support. This work is supported by a Grant-in-Aid for Scientific Research ”Novel Quantum Phenomena in Transition Metal Oxides” from the Ministry of Education, Science, Sports and Culture of Japan, Special Coordination Funds of the Science and Technology Agency of Japan, the New Energy and Industrial Technology Development Organization (NEDO), the U. S. DOE, Office of Basic Energy Science and Division of Material Science. SSRL is operated by the U. S. DOE, Office of Basic Energy Sciences, Division of Chemical Sciences.


text:
# ”Young” soap films

## 1 Introduction

Can a given surfactant produce a strong foam? there are some qualitative rules. We know, for instance, that a strongly insoluble surfactant cannot cover the surface of a rapidly growing bubble, while a soluble surfactant may succeed -by diffusion from the bulk solution. We would like to make this more precise- to characterise a surfactant solution by a few control parameters, which tell us what are its activities.

One of the methods for generating soap films is provided by a mechanical egg beater. Or, more scientifically, by pulling rapidly a metallic frame out of a surfactant solution (fig. ) as described in the book by Mysels et al .

We might first think of pulling the film at constant vertical speed. This is a version of a classical Landau Levich problem . But the quality of the surfactant plays only a minor role in this process : the surfactant is pulled up at (essentially) the same speed than the solution, and no fresh surface is generated.

A more relevant situation was invented by Lucassen . Here, the film is pulled abruptly, and its height $`h`$ is comparable to (or larger than) the horizontal span of the reservoir. Thus we must create fresh surface by pulling the frame, and migration of surfactant from the water phase to the surface is dominant.

This situation corresponds to what Lucassen  called a ”young film”. We discussed some aspects of the young films in a note . But we missed some important points: a) we chose as our central variable the local surface tension. It turns out that it is more illuminating to study first the concentration $`c(z)`$ at all levels in the film b) we did not discuss the dynamical features which tell us if the young film can indeed be made -or not c) some serious misprints occurred in .

Here, we return to the static structure of the young films: in section 2 we deal with concentrations $`c`$ below the vertical micelle concentration $`(c^{})`$. We find two characteristic heights: both are relevant for a discussion of film stability. In section 3 we extend this to $`c>c^{}.`$In section 4 we consider the major time constants involved: the young film must achieve its protective surface before falling.

## 2 Profile of a young film $`(c<c^{})`$

### 2.1 Basic equations

The film is drawn vertically (along $`z`$) and has a certain width $`\delta (z)`$ (fig. 1b). We assume that, in each interval $`dz,`$ there is a rapid equilibration between solution (concentration $`c(z)`$) and surface (with monolayer concentration $`\mathrm{\Gamma }(z)`$). The total number of surfactants was originally $`c_B\delta `$ (where $`c_B`$ is the bulk concentration in the reservoir), and remains the same. Thus we must have:

$$c_B\delta =\delta c+2\mathrm{\Gamma }$$
(1)
We assume that equilibrium has been achieved, and this may be written in the differential form:

$$\mathrm{\Gamma }^1d\mathrm{\Pi }=d\mu (c)$$
(2)
where $`\mathrm{\Pi }=\gamma _0\gamma `$ is the Langmuir pressure, describing the shift of the surface tension from $`\gamma _0`$ (for pure water) to $`\gamma `$ (at concentration $`c`$) and $`\mu (c)`$ is the chemical potential of the surfactant in solution.

These equations must be supplemented by a requirement of mechanical equilibrium: the weight of the film must be balanced by Marangoni forces:

$$\frac{2d\mathrm{\Pi }}{dz}=\rho g\delta $$
(3)
($`\rho =`$ water density, $`g=`$ gravitational acceleration).

### 2.2 Concentration profile

Combining eqs (1, 2, 3) we arrive at:

$$\frac{2d\mathrm{\Pi }}{dz}2\mathrm{\Pi }\frac{d\mu }{dz}\frac{dc}{dz}=\frac{\rho g2\mathrm{\Gamma }}{c_Bc}$$
(4)
In this section, we focus our attention on the dilute case $`c<c^{}`$. Then:

$$\mu =kT\mathrm{}nc+\text{cons}\mathrm{tan}\text{t}$$
(5)
and eq. 4 becomes:

$$\frac{dz}{\lambda }=\frac{dc}{c_B}\left(\frac{c_Bc}{c}\right)$$
(6)
where we have introduced a characteristic length $`\lambda `$ such that:

$$\rho g\lambda =c_BkT$$
(7)
$`\lambda `$ measures the osmotic pressure of the surfactant in terms of hydrostatic heights. Note that our definition of $`\lambda `$ differs from ref. . Eq. (6) integrates to:

$$\frac{z}{\lambda }=\frac{cc_B}{c_B}+\mathrm{}n\left(\frac{c_B}{c}\right)$$
(8)
This concentration profile is shown on fig. 2a. Two interesting limits are:

$$\frac{c_Bc}{c_B}=\left(\frac{2z}{\lambda }\right)^{1/2}(z<\lambda )$$
(9)
$$\frac{c}{c_B}=e^{z/\lambda }(z>\lambda )$$
(10)
### 2.3 Thickness profile

We can translate $`c(z)`$ into a thickness profile $`\delta (z)`$ using eq. (1). In most of the region $`c<c^{}`$, the monolayer concentration $`\mathrm{\Gamma }`$ is nearly constant (adsorption has taken place abruptly at very low $`c^{}s`$). Then we have:

$$\delta =\frac{2\mathrm{\Gamma }_B}{c_Bc}$$
(11)
where $`\mathrm{\Gamma }_B`$ is $`\mathrm{\Gamma }(c=c_B).`$ The thickness $`\delta (z)`$ described by eq. (11) is a decreasing function of $`c`$. When $`z>\lambda `$, we have $`c<<c_B,`$ and $`\delta `$ reaches a simple limit $`\delta =2\mathrm{}`$, where:

$$\mathrm{}=\mathrm{\Gamma }_B/c_B$$
(12)
The length $`\mathrm{}`$ plays an important role in foaming problems: $`\mathrm{}`$ is the minimal thickness of solution required to transfer surfactant from bulk to surface, achieving the required surface concentration $`\mathrm{\Gamma }_B.`$

### 2.4 Vertical variations of the surface tension $`\gamma `$

We know $`c(z)`$, and we know (from classical plots for usual surfactants) the surface tension $`\gamma (c)`$: thus we may construct $`\gamma (z)`$=$`\gamma _0\mathrm{\Pi }(z).`$ The general aspect of this curve is shown on fig. 2b.

For $`z>\lambda `$ ($`c<<c_B`$), we may write:

$$\frac{d\gamma }{dz}=\frac{d\mathrm{\Pi }}{dz}=\frac{\rho g\mathrm{\Gamma }}{c_Bc}\frac{\rho g\mathrm{\Gamma }_B}{c_B}=\text{constant}$$
(13)
Thus $`\gamma (z)`$ increases linearly with $`z`$, up to a certain height ($`h_m`$), where we return to a bare surface ($`\mathrm{\Gamma }0,`$ $`\gamma =\gamma _0`$). This corresponds to:

$$h_m\frac{d\gamma }{dz}\gamma _0\gamma _B\mathrm{\Pi }_B$$

or:

$$h_m=\lambda \frac{\mathrm{\Pi }_B}{kT\mathrm{\Gamma }_B}$$
(14)
Note that $`h_m`$ is larger than $`\lambda `$,because $`\mathrm{\Pi }_B`$ si much larger than an ideal gas pressure. The height $`h_m`$ is an absolute limit of stability for our young film. (But of course, other instabilities may occur in the regime of high $`z`$, low $`\delta `$).

## 3 Extension to higher concentrations $`c_B>c^{}`$

### 3.1 Structure of the chemical potential at $`c>>c^{}`$

At concentrations $`c>c^{}`$, most of the surfactant is in a micellar form (concentration $`c_M`$). Only a small fraction ($`c_1`$) is present in a monomer form. Let us call $`N`$ the number of surfactants per micelle. Then a rough but convenient description (with $`N`$ fixed) of the micelle/monomer equilibrium, may be written in the form:

$$c_M=\frac{c_1^N}{c^{N1}}$$
(15)
or for the total concentration:

$$c=c_1+c_M=c_1\left\{1+\left(\frac{c_1}{c^{}}\right)^{N1}\right\}$$
(16)
The chemical potential is still:

$$\mu =kT\mathrm{}nc_1+\text{cons}\mathrm{tan}\text{t}$$
(17)
Let us consider only the limit $`c>>c^1`$. Then:

$$\frac{c}{c^1}\frac{c_M}{c^{}}=\left(\frac{c_1}{c^{}}\right)^N$$
(18)
and:

$$\mu =kT\left\{\mathrm{}nc^{}+\frac{1}{N}\mathrm{}n\frac{c}{c^{}}\right\}$$
(19)
$$\frac{d\mu }{dc}=\frac{kT}{Nc}$$
(20)
### 3.2 Concentration profile

Returning to eq. (4), we see from eq. (20) that for $`c>>c^{},`$ we should replace $`\rho g`$ by $`N\rho g.`$ Or, equivalently, the characteristic length in the high concentration regime is:

$$\stackrel{}{\lambda }=\frac{\lambda }{N}$$
(21)
This is much smaller than $`\lambda `$ since $`N`$ is large ($``$ 60).

The resulting plot of $`c(z)`$ is shown on fig. 3a. We start at $`z=0`$ from a concentration $`c_B>>c^{}.`$Then upon increasing $`z`$, the concentration decreases, and when $`c<<c_B`$, we have simply:

$$c(z)=c_B\mathrm{exp}\left(\frac{Nz}{\lambda }\right)(z<z^{})$$
(22)
This holds up to a height $`z^{}`$ where $`c=c^{}:`$

$$z^{}=\frac{\lambda }{N}\mathrm{}n\frac{c_B}{c^{}}$$
(23)
At $`z>z^{}`$, we return to the low concentration regime (eq. 6), and ultimately, for $`z>\lambda `$, we have:

$$c(z)c^{}\mathrm{exp}\left(\frac{z}{\lambda }\right)$$
(24)
### 3.3 Surface tensions

In all the region $`c>c^{},`$ or $`z<z^{}`$, the surface tension keeps its minimum value $`\gamma ^..`$Above $`z<z^{},`$ we return to the cons$`\mathrm{tan}`$t slope regime of eq. (13). The maximum height allowable $`h_m`$ is still given by eq. (14) as shown on fig. (3b).

Note the inequalities:

$$h_m>\lambda >z^{}$$
(25)
## 4 Discussion

### 4.1 Summary of static predictions:

a) in a young film, we expect a concentration $`c(z)`$ which decays exponentially at large $`z`$, with a characteristic length $`\lambda `$ (eq. 7). $`\lambda `$ is typically a few meters. This statement holds for bulk concentrations $`c_B`$ which may be lower or higher than $`c^{}`$.

b) the film thickness at high $`z`$ is simply $`\delta =2\mathrm{}`$, where $`\mathrm{}=\mathrm{\Gamma }_B/c_B`$ -$`(c<c^{})`$ or $`\mathrm{}=\mathrm{\Gamma }_B/c^{}`$ for $`c>c^{}.`$

c) the film is bare and completely unstable at heights $`z>h_m`$, where $`h_m`$ is given by eq. (14), and is larger than $`\lambda .`$

d) on the whole, going to bulk concentrations, $`c_B>c^{}`$ does not produce great alterations: most static film properties should be close to what they are for $`c_B=c^{}.`$

### 4.2 Characteristic times

1) A water film with no surfactant will collapse under its own weight in a time $`t_c`$ governed by gravitational flow with the acceleration $`g:`$

$$t_c\sqrt{\frac{h}{g}}$$
(26)
where $`h`$ is the height of the film. For $`h=10`$ cm, $`t_c0.1`$ sec.

2) What is the time required for the surfactant to reach the surface? We consider here first the (relatively simple) case where $`c_B<c^{}`$: all the transport is via monomers. We also assume that there is no barrier at the surface opposing the surfactant adsorption. Then, we must transfer the surfactant by diffusion (coefficient $`D`$) over a distance $`\mathrm{}`$ defined in eq. (12). The diffusion time is roughly:

$$t_d=\frac{\mathrm{}^2}{D}=\frac{\mathrm{\Gamma }_B^2}{c_B^2D}$$
(27)
The condition of formation of a young film at $`c<c^{}`$ is essentially that $`t_d`$ be shorter than $`tc`$. To get small $`t_d^{}s`$, we must use concentrations $`c_B`$ which are as high as possible (but remain below $`c^{}`$).

3) The extension of these ideas to $`c_B>c^{}`$ is quite delicate. Micelles provide a big reservoir of surfactant, but the delivery is slow.

a) micelles can emit monomers which are then absorbed by the fresh surface: this process is fast (10-6 seconds) but inefficient, because each micelle cannot give more than 2 or 3 monomers (beyond that, the micellar free energy rises fast).

b) micelles can split in two, and the corresponding characteristic time $`\tau _{split}`$ is relatively long (10-3 sec). After splitting, the parts end up as monomers, and can feed the surface. This process has been documented by Shah .

c) the micelles may unfold by direct contact with the surface. The energy barriers for processes (b) and (c) are comparable -maybe the barrier for (c) is slightly lower.

If, for simplicity, we focus on process (b), the conclusion is that if we have $`t_{split}<t_c`$, the young film should be realisable. The opposite limit ($`t_{split}>t_c`$) would occur only for very long surfactant chains.

### 4.3 Conclusions

a) To achieve young films of height $`h`$, we need $`h_m>h`$, where $`h_m`$ is given by eq. (14). As shown by this equation, high $`h_m`$ corresponds to high surface pressures (or $`\gamma <<\gamma _0`$) in the present solution. To achieve this, our surfactants must have an aliphatic tail which is not too short.

b) The surfactant must migrate fast to the water/air interface. For $`cc^{}`$, this implies that $`\mathrm{}^{}=\mathrm{\Gamma }^{}/c^{}`$ must be relatively small. The critical micelle concentration $`c^{}`$ should not be too small. The surfactant tails should not be too long.

Thus a good foaming agent results from a compromise: this does correspond to the empirical rule stating that the HLB should be in the range 7 to 9.

Acknowledgments: we have benefited from useful advice by D. Langevin and M. Schott.


text:
# Topological Description of the Aging Dynamics in Simple Glasses

\[

## Abstract

We numerically investigate the aging dynamics of a monatomic Lennard-Jones glass, focusing on the topology of the potential energy landscape which, to this aim, has been partitioned in basins of attraction of stationary points (saddles and minima). The analysis of the stationary points visited during the aging dynamics shows the existence of two distinct regimes: i) at short times, $`t<t_c`$, the system visits basins of saddles whose energies and orders decrease with $`t`$; ii) at long times, $`t>t_c`$, the system mainly lies in basins pertaining to minima of slowly decreasing energy. The dynamics for $`t>t_c`$ can be represented by a simple random walk on a network of minima with a jump probability proportional to the inverse of the waiting time.

\]

The very slow approach to equilibrium of glassy systems is a very interesting and intriguing physical process that, in the last few years, has received considerable attention from analytical, numerical and experimental point of view . When a glassy system is suddenly brought in the glassy state from an high temperature (or low density) equilibrium state, it starts to explore the phase space with an off-equilibrium dynamics which slows down as the system “ages”. During this aging regime one time physical quantities, that are time independent on average at equilibrium, slowly approach their equilibrium values, while two time quantities, like the correlation functions, depend on both times due to the lack of time translation invariance. The study of the off-equilibrium dynamics of glassy systems, and in particular of spin-glasses, has allowed to obtain important informations on the phase-space properties of the systems itself (behavior of the order parameter in the low temperature equilibrium phase ) and to formulate off-equilibrium relations generalizing the usual equilibrium ones (generalized fluctuation-dissipation relation ) with the introduction of an effective temperature . The investigation of the off-equilibrium properties seems to indicate the correctness of the conjecture stating the similarity of structural glasses with some spin glass model (spin-glasses with one step replica symmetry breaking) . The way a glassy system goes towards equilibrium is then not only an interesting problem by itself but it is also of great relevance in view of a deeper understanding of the nature of glass transition. In this respect, many efforts have been devoted to the understanding of the role played by the potential energy landscape in the equilibrium dynamics of supercooled liquids . The trajectory of the representative point in the $`3N`$ configuration space can be mapped into a sequence of locally stable points (the so called Inherent Structures, IS ), that are the local minima of the total potential energy $`V`$: to each instantaneous configuration during the dynamical evolution of the system one can associate an IS by a steepest descent path in the $`V`$ surface. The properties of the IS has been found to be very useful to clarify many features in the dynamics and the thermodynamics of supercooled liquids in  and off  equilibrium. An extension of such an analysis including all stationary points of $`V`$ (minima and saddles)  has been recently proposed to form a better representation of the system during the equilibrium dynamics. This approach gives a novel point of view on the equilibrium dynamics of supercooled liquids based on a detailed topological analysis of the potential energy landscape. In particular it emerges  that the order of the sampled saddles during the equilibrium dynamics is a well defined increasing function of temperature and it appears to vanish at the so-called mode-coupling  transition temperature $`T_{MCT}`$, which then marks the crossover from a dynamics between basins of saddles ($`T>T_{MCT}`$) to a dynamics between basins of minima ($`T<T_{MCT}`$). Similar results have been found  in the framework of Instantaneous Normal Modes analysis .

In this Letter we numerically investigate the off-equilibrium dynamics of a simple model glass, by mapping the true dynamics in the configuration space on to the stationary points of $`V`$. In our glass the particles interact via a Lennard-Jones potential modified in such a way as to avoid crystallization . The system is brought in an off-equilibrium initial condition by a sudden isothermal density jump (crunch) across the liquid-glass transition line. We find a rather sharp cross-over between two dynamical regimes: i) a short times regime where the system point travels between basins of saddles of decreasing energy and order and ii) a long time regime where the system point moves between basins of saddles of zero order (minima) and decreasing energy, crossing basins of saddles of low order. The existence of two well separated dynamical regimes shows up also in the Mean Square Displacement (MSD). Here a logarithmic behavior is found at long times suggesting the picture of a system point performing a random walk between basins of minima with a jump probability per unit time inversely proportional to time elapsed since crunch. This picture is quantitatively supported by the very good agreement found between the MSD measured along the MD trajectories and the MSD reconstructed via a random walk model.

The system under study is a monatomic Modified-Lennard-Jones (MLJ) that is the usual 6-12 Lennard-Jones pair interaction (truncated and shifted at $`R_c=2.6`$, standard Lennard-Jones units are used hereafter) plus a small many body term which inhibits crystallization with negligible corrections to the equation of state of a monatomic Lennard-Jones system (see  for details) . The system is composed of $`N=256`$ particles enclosed in a cubic box with periodic boundary conditions. The initial liquid equilibrium configurations ($`\rho _0=0.95`$, $`T_0=0.5`$) are prepared by standard isothermal MD. At time $`t=0`$ the density is suddenly increased to the value $`\rho _1=1.24`$, where the glass transition temperature ($`T_g(\rho _1)1.4`$ ) is well above the simulation temperature $`T_0`$. The subsequent evolution of the system is generated up to time $`t=310^4`$ by means of standard microcanonical MD with time step $`\tau =0.01`$. For each MD trajectory we quenched a set of logarithmically equally spaced configurations to the corresponding IS and saddle following the steepest descent path respectively on the potential energy surface $`V`$ in the former case and on the auxiliary potential $`W=\frac{1}{2}|\stackrel{}{}V|^2`$  in the latter.

In Fig.1 we report in a logarithmic scale the time evolution of: a) the instantaneous configuration potential energy per particle $`e_i`$ subtracted of the constant term $`\frac{3}{2}T_0`$ (open circles), b) the saddle’s energy $`e_s`$ (full triangles), c) the IS energy $`e_m`$ (open diamonds). The remarkable coincidence between $`e_s`$ and $`e_i\frac{3}{2}T_0`$ suggests to think of the instantaneous energy $`e_i`$ as the sum of two contributions: the saddle’s energy (slowly approaching equilibrium) and the energy of nearly $`N`$ fast harmonic equilibrated degrees of freedom. It is evident the presence of a rather sharp transition in dynamical behavior of $`e_i`$ and $`e_s`$ at time $`t_c20`$. In particular the time evolution of $`e_i`$ for $`t>t_c`$ seems to be driven by the inherent dynamics between minima as suggested by the underlying time evolution of $`e_m`$.

In Fig.2 the mean saddle’s order $`n_s`$ (number of negative eigenvalues of the potential energy Hessian computed on the saddle configurations) is reported as a function of time. The dynamic transition observed at $`t_c`$ in the energy time evolution in Fig.1 corresponds here in the mean order $`n_s`$ falling below the value $`n_s=1`$ or, in other words, for times longer than $`t_c`$ the system point starts to explore extensively the basins of saddles of order zero (minima). The inset in Fig.2 shows the order $`n_s`$ computed along a single MD trajectory. One can see that for times $`t>t_c`$ the systems visits basins of minima passing through basins of low order saddles.

In a recent study of the MLJ landscape probed in equilibrium  some of us found that the mean elevation energy of a saddle of order $`n_s`$ from the underlying minimum is proportional to $`n_s`$ (this is at variance with the result reported in , where a linear relation is observed between $`n_s`$ and the absolute energy of the saddles). This is clearly a property of the energy landscape itself and one expects that the same properties are also valid in the off-equilibrium dynamics.

In Fig.3 where we report $`e_s`$ together with $`e_m+n_s\mathrm{\Delta }E/N`$ with $`\mathrm{\Delta }E`$ a fit parameter. This result confirms the picture of an energy landscape where the saddles of order $`n_s+1`$ lay above the saddles of order $`n_s`$ of a constant quantity $`\mathrm{\Delta }E`$. From the data we obtain the value $`\mathrm{\Delta }E=10`$. The value of $`\mathrm{\Delta }E`$ could also be obtained from an equilibrium analysis, as in  where a value of $`\mathrm{\Delta }E=3.6`$ is found for the case $`\rho =1`$. We have not carried out such an equilibrium analysis for the density $`\rho =1.24`$ used here.

Another interesting quantity in describing the dynamical processes is the mean square displacement of the configuration at time $`t`$ from that at time $`t_w`$ after crunch:

$$r^2(t,t_w)=\frac{1}{N}\underset{i}{}|𝐫_i(t)𝐫_i(t_w)|^2,$$
(1)
where $`\mathrm{}`$ is an average over initial configurations. In Fig.4 we report $`r^2(t,t_w)`$ as a function of time $`t`$ for four selected $`t_w`$ values: $`t_w=1,5,20,200`$. Also in this case $`t_c20`$ marks a transition in the dynamical behavior. In particular the long time regime shows a well defined logarithmic dependence on time. This suggests the picture of the diffusion process proceeding through elementary uncorrelated events occurring with a probability per unit time $`p(t)1/t`$ . As previously observed, in the long time regime the system point travels between basins of minima. It is therefore natural to identify the elementary diffusion events as jumps between adjacent basins of minima. In order to give further support to this conjecture we calculated the cumulative number of jumps $`N_{jumps}(t)`$ defined as the number of jump events in which the constantly monitored energy of the IS makes a sharp transition.

The quantity $`N_{jumps}`$ is plotted in Fig.5(a). For times $`t>t_c`$ $`N_{jumps}`$ displays a nice logarithmic behavior that is what one would expect if the probability per unit time of a jump event is $`p(t)=c/t`$. In other words the time needed to escape from an IS reached by the system after a time $`t_w`$ from crunch is proportional to time $`t_w`$ itself (a similar behavior is found in the off-equilibrium dynamics of a ferromagnetic Ising system where the persistence time of a bubble in which spins have a similar value is proportional to time $`t_w`$ ). In Fig.5(a) the solid line represents the fitting function $`c\mathrm{log}(t)`$ with $`c=1.0`$. In Fig.5(b) we report the mean square distance $`d^2`$ between adjacent minima probed by the system during the aging dynamics:

$$d^2(t)=\frac{1}{N}\underset{i}{}|𝐫_i^{IS}(t^+)𝐫_i^{IS}(t^{})|^2,$$
(2)
where $`t`$ is the time of the jump between two different basins of minima and $`𝐫_i^{IS}(t^{})`$, $`𝐫_i^{IS}(t^+)`$ are the coordinates of IS respectively before and after time $`t`$. The quantity reported in figure approaches the long time constant value $`d_0^20.019`$. For uncorrelated jumps, the mean square displacement would be simply given by

$$r^2(t,t_w)=d_0^2c\mathrm{log}(t/t_w).$$
(3)
This is in a remarkably good agreement with the MD data as shown in Fig.4, where we draw straight lines with slope $`d_0^2c=0.019`$ (as obtained from the values of $`d_0^2`$ and $`c`$ derived from Fig.5) and intercept chosen as to best fit the data for $`t>20`$.

In conclusion we have numerically studied the off-equilibrium dynamics of a monatomic Modified-Lennard-Jones system by sudden isothermal density jumps from the liquid to the glassy phase. Analyzing the time evolution of the energy and the mean square displacement we found two well defined dynamical regimes, which can be related to the different ways the system explores the potential energy landscape: i) at short times the system point moves downward passing through basins of saddles with decreasing energy and order, and ii) at longer time the relevant dynamical processes are jumps between basins of minima with slower decreasing energy. In the whole time regime the instantaneous potential energy $`e_i`$ can be thought as due to a sum of two contributions: the energy of the underlying saddle’s configurations $`e_s`$ plus a nearly harmonic term of the fast degrees of freedom $`3/2T`$. Moreover the linear energy gap between saddles and underlying minima $`e_se_m`$ as a function of saddles order $`n_s`$ supports the picture of a simple energy landscape with only a single energy barrier parameter $`\mathrm{\Delta }E`$ which represents the energy gap between a saddle of order $`n_s`$ and a saddle of order $`n_s+1`$. The relevant processes in the long time aging regime are well represented by a simple random walk model with elementary events corresponding to jumps between minima with a probability per unit time proportional to the inverse of time elapsed since the system has been brought out of equilibrium.


text:
# On the spin–up of neutron stars to millisecond pulsars in long–period binaries

## 1. Introduction

Ever since the detection of the first millisecond pulsar (Backer et al. 1982) it has been clear that a ms-pulsar is a neutron star (NS) which has been spun up by accretion in a close binary system (for a review see Bhattacharya & van den Heuvel 1991). Because the amount of mass needed to spin up a NS to ms periods is rather small, i.e. $`0.05\mathrm{M}_{}<\mathrm{\Delta }M_{\mathrm{NS}}<0.1\mathrm{M}_{}`$ (Burderi et al. 1999) and the mass available from the donor is typically many times larger, i.e. typically of order $`1\mathrm{M}_{}`$, it appeared as if there was no problem spinning up a NS to such short periods. Yet in recent years, evidence has accumulated which shows that NS are inefficient accretors. Thorsett & Chakrabarty (1999) noted that in all cases where the mass of a NS could be measured with precision it turned out to be in the very narrow interval $`1.3\mathrm{M}_{}<M_{\mathrm{NS}}<1.4\mathrm{M}_{}`$. Even if one makes the extreme assumption that all NSs are born with exactly the same mass, this would imply that none of these NSs could have accreted more than $`0.1\mathrm{M}_{}`$, despite the fact that some of these NSs have been spun up by accretion and the mass available from the donor star was probably much larger than $`0.1\mathrm{M}_{}`$. In another context, King & Ritter (1999) have presented further evidence that NSs can be very inefficient accretors. The main reason for this low accretion efficiency is the frequent disparity between the mass transfer rate provided by the donor star and the maximum rate at which the NS can accrete. The latter is approximately the Eddington rate $`\dot{M}_\mathrm{E}1.5\times 10^5\mathrm{M}_{}\mathrm{yr}^1`$, whereas, depending on circumstances, the former can be larger by many orders of magnitude.

One of the obvious places in which NSs can be spun up to ms–pulsars are wide binaries in which a NS accretes from a giant. Tauris & Savonije (1999) have recently devoted a detailed study to this formation channel. One of their predictions is that, in principle, ms–pulsars should form in binaries with final orbital periods of up to $`P1000^\mathrm{d}`$, even when taking into account that the mass transfer rate in long–period systems ($`P>20^\mathrm{d}`$) exceeds the Eddington accretion rate. In addition, their calculations predict a correlation between the final orbital period and the mass of the ms–pulsar. Unfortunately, this cannot be tested at present. But their first prediction can. And here observations show, contrary to their prediction, that there are no ms–pulsars in wide binaries with a circular orbit and a period $`P>200^\mathrm{d}`$ (see e.g. table 1 in Taam, King & Ritter 2000). In this context, Li & Wang (1998) have already noted that the accretion efficiency of the NS in such a binary is far lower than has been assumed by Tauris & Savonije (1999) if one takes into account the fact that in all these wide binaries the accretion disc around the NS is thermally unstable and accretion is transient (King et al. 1997). This is because a) during an outburst the mass flow rate through the disc is much larger than the mass transfer rate from the donor, which, in turn, is already super–Eddington if $`P>20^\mathrm{d}`$, and b) because a NS cannot accrete in excess of the Eddington rate. In fact, Li & Wang (1998) found that if the duty cycle of the disc instability limit cycle is $`d<10^2`$ then in systems with a final orbital period $`P_\mathrm{f}>100^\mathrm{d}`$ the total mass accreted by the NS is $`\mathrm{\Delta }M_{\mathrm{NS}}<`$ few $`10^2\mathrm{M}_{}`$, too little to spin up a NS to ms periods.

In this paper we elaborate on the idea of Li & Wang (1998) and derive explicit expressions for the duty cycle and the accretion efficiency of the NS in terms of the parameters of the binary system and the outburst cycle. Our calculations show that if the duration of the quiescent phase between two subsequent outbursts exceeds a few decades, the accretion efficiency becomes so small that in binaries with a final orbital period $`P_\mathrm{f}>200^\mathrm{d}`$, the NS cannot accrete enough mass to become a ms pulsar.

## 2. Basic model ingredients

The binary systems for which we wish to decribe the accretion processes are wide, long–period binaries with orbital periods of typically $`P>20^\mathrm{d}`$ and in which a neutron star primary of mass $`M_1`$ accretes from a giant secondary, of mass $`M_2`$. The basic model ingredients we use can be summarized as follows:

### 2.1. Mass transfer

We assume that mass transfer from the secondary is thermally (and adiabatically) stable. In this case, the mass transfer rate can be computed from a simple analytical approximation (e.g. Ritter 1999). Following King et al. (1997) and using the same model parameters we can write for the mass transfer rate

$$\dot{M}_2=7.3\times 10^{10}\mathrm{M}_{}\mathrm{yr}^1(\zeta _e\zeta _R)^1m_2^{1.7425}p^{0.9281},$$
(1)
where $`m_2=M_2/\mathrm{M}_{},p=P(\mathrm{d})`$ is the orbital period in days, $`\zeta _e`$ the thermal equilibrium mass–radius exponent of the donor star, and $`\zeta _R`$ the mass–radius exponent of the donor’s critical Roche radius. For (thermal) stability we require $`\zeta _e\zeta _R>0`$ (e.g. Ritter 1988).

### 2.2. Disc irradiation

Accretion on to a compact star can result in significant external irradiation of the accretion disc by the central accretion source (e.g. van Paradijs 1996). Here we follow King & Ritter (1998) and assume that the effective temperature of the disc is kept above the hydrogen ionization temperature $`T_\mathrm{H}6500K`$ by irradiation from the central accretion source out to a radius

$$R_\mathrm{h}=(B_1\dot{M}_c)^{1/2},$$
(2)
where $`\dot{M}_c`$ is the central accretion rate, and, for a neutron star accretor, $`B_1=3.9\times 10^5`$ cgs (King & Ritter 1998). Furthermore, we shall assume that the neutron star cannot accrete in excess of the Eddington accretion rate $`\dot{M}_\mathrm{E}1.5\times 10^8\mathrm{M}_{}\mathrm{yr}^1`$, i.e. that $`\dot{M}_c<\dot{M}_\mathrm{E}`$. This means that the maximum radius out to which the disc can be kept hot by irradiation, i.e. at $`T_{\mathrm{eff}}>T_\mathrm{H}`$, is

$$R_{\mathrm{crit}}7.1\mathrm{R}_{}(b_1\dot{m}_\mathrm{E})^{1/2}.$$
(3)
Here $`\dot{m}_\mathrm{E}=\dot{M}_\mathrm{E}/10^8\mathrm{M}_{}\mathrm{yr}^1`$ and $`b_1=B_1/3.9\times 10^5`$ cgs. At the same time our assumption that $`\dot{M}_c<\dot{M}_\mathrm{E}`$ means that whenever the mass flow rate through the disc is $`\dot{M}_\mathrm{d}>\dot{M}_\mathrm{E}`$, the excess matter must be ejected from the disc.

### 2.3. The $`\alpha `$–disc model

For describing the properties of the accretion disc around the neutron star we use the Shakura & Sunyaev (1973) $`\alpha `$–disc model. Accordingly, if the vertical temperature stratification of the disc is dominated by viscous heating, the disc viscosity is (e.g. Frank, King & Raine 1992)

$$\nu _{\mathrm{visc}}=2.7\times 10^{15}\mathrm{cm}^2\mathrm{s}^1\alpha ^{4/5}(\dot{m}_\mathrm{d})^{3/10}m_1^{1/4}r^{3/4},$$
(4)
where $`\dot{M}_\mathrm{d}`$ is the mass flow rate through the disc, $`R`$ the disc radius, $`\dot{m}_\mathrm{d}=\dot{M}_\mathrm{d}/10^8\mathrm{M}_{}\mathrm{yr}^1`$, and $`r=R/\mathrm{R}_{}`$. $`\alpha `$ is the viscosity parameter, here assumed to be constant with radius.

(4) holds only to the extent that the vertical temperature stratification in the disc is dominated by viscous heating, hence its dependence on $`\dot{M}_\mathrm{d}`$. If, on the other hand, the disc is strongly heated by external irradiation, the vertical temperature profile becomes nearly isothermal at the irradiation temperature $`T_{\mathrm{irr}}`$. Taking this into account in the $`\alpha `$–viscosity ansatz we find (see also King 1998)

$$\begin{array}{cc}\nu _{\mathrm{visc}}\hfill & =\alpha \frac{T_\mathrm{H}}{\mu }(B_1\dot{M}_c)^{1/4}(GM_1)^{1/2}R\hfill \\ & =3.8\times 10^{15}\mathrm{cm}^2\mathrm{s}^1\alpha \left(\frac{\mu }{0.6}\right)^1\left(\frac{T_H}{6500K}\right)(b_1\dot{m}_c)^{1/4}m_1^{1/4}r.\hfill \end{array}$$
(5)
Here $``$ is the gas constant and $`\mu `$ the mean molecular weight.

The maximum suface density which a disc can reach in the cool state, i.e. before it becomes thermally unstable and must turn into the hot state, is conventionally denoted by $`\mathrm{\Sigma }_B`$. Taking $`\mathrm{\Sigma }_B`$ from computations of Ludwig, Meyer–Hofmeister & Ritter (1995), the mass accumulated before an outburst sets in can be written as

$$M_d7\times 10^{10}\mathrm{M}_{}\alpha _c^{4/5}m_1^{0.37}r_\mathrm{d}^{3.10}.$$
(6)
Here $`r_\mathrm{d}=R_\mathrm{d}/\mathrm{R}_{}`$, where $`R_\mathrm{d}`$ is the outer radius of those parts of the disc which are involved in the outburst, $`\alpha _c`$ the $`\alpha `$–parameter in the cool state, and $`<1`$ is the filling factor defined by this equation. $``$ takes into acount that a the onset of an outburst $`\mathrm{\Sigma }=\mathrm{\Sigma }_B`$ at only one point in the disc and that $`\mathrm{\Sigma }<\mathrm{\Sigma }_B`$ elsewhere.

### 2.4. Assumptions about the disc instability limit cycle

We shall assume that the disc instability limit cycle in these wide, long–period binary systems with a neutron star accretor can be described as follows:

During an outburst the disc is kept in the hot state by irradiation out to the radius $`R_\mathrm{h}(\dot{M}_c)`$ given by (2). Because $`\dot{M}_c<\dot{M}_\mathrm{E},R_\mathrm{h}(\dot{M}_c)<R_\mathrm{h}(\dot{M}_\mathrm{E})`$. If the physical outer radius of the disc $`R_{\mathrm{disc}}>R_\mathrm{h}(\dot{M}_c)`$, only the inner part with radius $`R_\mathrm{d}=R_\mathrm{h}(\dot{M}_c)`$ is hot, whereas the outer part with $`R_\mathrm{h}(\dot{M}_c)<R<R_{\mathrm{disc}}`$ remains in the cold state (e.g. King & Ritter 1998). Hereafter we shall call $`R_\mathrm{d}`$ the radius of the “active” disc. If, on the other hand, $`R_\mathrm{h}(\dot{M}_c)>R_{\mathrm{disc}},R_\mathrm{d}=R_{\mathrm{disc}}`$. If the mass flow rate through the outer parts of the active disc is $`\dot{M}_\mathrm{d}(R_\mathrm{d})>\dot{M}_\mathrm{E}`$, then as long as $`\dot{M}_\mathrm{d}(R_\mathrm{d})>\dot{M}_\mathrm{E},\dot{M}_c=\dot{M}_\mathrm{E},`$ and $`R_\mathrm{h}(\dot{M}_c)=R_\mathrm{h}(\dot{M}_\mathrm{E})`$ = const.

During an outburst the active part of the disc (out to $`R_\mathrm{d}`$) is essentially emptied, i.e. the mass remaining in the active part of the disc immediately after an outburst is insignificant compared to the mass $`M_\mathrm{d}`$ given by (6) immediately before the onset of an outburst.

During quiescence the mass of the active part of disc which was accreted/ ejected during the previous outburst is replenished at the mass transfer rate $`\dot{M}_2`$. This is tantamount to assuming that if an inactive outer part of the disc exists, i.e if $`R_\mathrm{h}(\dot{M}_c)<R_{\mathrm{disc}}`$, this outer part is stationary, i.e. that $`\dot{M}_\mathrm{d}(R)=\dot{M}_2`$ for all $`R_\mathrm{h}(\dot{M}_c)<R<R_{\mathrm{disc}}`$.

We assume that the hot part of the disc (during an outburst) is characterized by a viscosity parameter $`\alpha _\mathrm{h}`$, whereas in the cool state (i.e. in the active part during quiescence and in the inactive part) the viscosity parameter is $`\alpha _c<\alpha _\mathrm{h}`$.

## 3. Disc accretion in long–period neutron star binaries

Here we are mainly interested in systems which end their nuclear time scale–driven mass transfer at very long orbital periods, i.e. $`P>200^\mathrm{d}`$. With the analytical solution (Ritter 1999) one can show that the corresponding initial periods at which stable mass transfer must have started are typically $`P_\mathrm{i}>20^\mathrm{d}`$, (the precise value depending on whether mass transfer is conservative or not, cf. Taam et al. 2000). It has been shown by King et al. (1997) and Li & Wang (1999) that at such long periods the mass transfer rate provided by nuclear evolution of the donor is always too small for stable disc accretion, i.e. the systems in question are transient, going through a disc instability limit cycle.

In addition, the orbital period at which $`R_{\mathrm{disc}}=R_{\mathrm{crit}}`$ (with $`R_{\mathrm{disc}}>R_{\mathrm{crit}}`$ for longer periods) is

$$P_{\mathrm{crit}}=12\mathrm{d}\left(\frac{R_{\mathrm{disc}}/R_{1,R}}{0.7}\right)^{3/2}m_1^{0.675}m_2^{0.175}(b_1\dot{m}_\mathrm{E})^{3/4},$$
(7)
where $`R_{1,R}`$ is the critical Roche radius of the primary. Thus, for practically all systems in question we also have $`R_{\mathrm{disc}}>R_\mathrm{h}(\dot{M}_c)`$.

As we shall argue below, during most of an outburst the mass flow rate through the outer parts of the active disc is highly super–Eddington in the sense that $`\dot{M}_\mathrm{d}(R_\mathrm{d})\dot{M}_\mathrm{E}`$. Therefore, at least initially, in these systems $`R_\mathrm{d}=R_\mathrm{h}(\dot{M}_\mathrm{E})`$ = const.

Now, if the viscosity at $`R=R_\mathrm{d}`$ does not depend on $`\dot{M}_\mathrm{d}`$, and $`R_\mathrm{d}`$ itself is constant, both the mass and the mass flow rate through the outer parts of the active disc decay exponentially on the viscous time scale

$$\tau _{\mathrm{visc}}=\frac{R_\mathrm{d}^2}{3\nu _{\mathrm{visc}}(R_\mathrm{d})}$$
(8)
at the outer radius of the active disc. If initially $`\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)\dot{M}_\mathrm{E}`$, we have $`R_\mathrm{d}=R_\mathrm{h}(\dot{M}_\mathrm{E})=\mathrm{const}.`$ as long as $`\dot{M}_\mathrm{d}>\dot{M}_\mathrm{E}`$ and the disc is essentially emptied during the super–Eddington phase which lasts for a time

$$t_{\mathrm{outb}}\tau _{\mathrm{visc}}(R_\mathrm{d})\mathrm{ln}\left[\frac{\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)}{\dot{M}_\mathrm{E}}\right].$$
(9)
This is the behaviour to be expected if the active part of the disc is dominated by irradiation, i.e. if $`\nu _{\mathrm{visc}}`$ is given by (5). If, on the other hand, viscous heating still dominates the vertical temperature structure, $`\nu _{\mathrm{visc}}`$ is given by (4). Using (4) in (8) we see that $`\dot{M}_\mathrm{d}M_d^{10/7}`$, i.e. that the mass flow rate and the mass in the active part of the disc do not decay exponentially but rather as

$$M_d(t)=M_d(0)\left[1+\frac{3}{7}\frac{t}{\tau _0}\right]^{3/7},$$
(10)
where

$$\tau _0=\frac{M_d(0)}{\dot{M}_\mathrm{d}(0)}$$
(11)
is the time scale on which the disc mass decays initially. If, initially, $`\dot{M}_\mathrm{d}\dot{M}_\mathrm{E}`$, the disc is essentially emptied during the super–Eddington phase which lasts in this case for the time

$$t_{\mathrm{outb}}\frac{3}{7}(\tau _{\dot{M}_\mathrm{E}}\tau _0),$$
(12)
where

$$\tau _{\dot{M}_\mathrm{E}}=\tau _0\left[\frac{\dot{M}_\mathrm{d}(0)}{\dot{M}_\mathrm{E}}\right]^{3/10}$$
(13)
is the time scale on which the disc mass decreases once $`\dot{M}_\mathrm{d}`$ reaches $`\dot{M}_\mathrm{E}`$.

The duration of the quiescent phase is given by the time it takes to replenish the matter which has been lost (accreted/ejected) from the active part of the disc during the previous outburst. Because we have assumed that the active disc is fed at the rate $`\dot{M}_2`$, irrespective of whether $`R_{\mathrm{disc}}`$ is larger or smaller than $`R_\mathrm{d}`$, the duration of quiescence is

$$t_\mathrm{q}\frac{M_d(0)}{\dot{M}_2}.$$
(14)
The total length of an outburst cycle is then

$$t_{\mathrm{cycle}}=t_{\mathrm{outb}}+t_\mathrm{q}$$
(15)
and the duty cycle of the outbursts

$$d=\frac{t_{\mathrm{outb}}}{t_{\mathrm{cycle}}}.$$
(16)
To decide the question of how much the neutron star can accrete, the quantity we need is the accretion efficiency $`\eta `$. Because the mass accreted over one cycle is $`\mathrm{\Delta }M_{\mathrm{accr}}t_{\mathrm{outb}}\dot{M}_\mathrm{E}`$, whereas the mass transferred over the same time is $`\mathrm{\Delta }M_{\mathrm{transf}}t_{\mathrm{cycle}}(\dot{M}_2)`$, we have

$$\eta \frac{t_{\mathrm{outb}}\dot{M}_\mathrm{E}}{t_{\mathrm{cycle}}(\dot{M}_2)}d\frac{\dot{M}_\mathrm{E}}{(\dot{M}_2)}.$$
(17)
Writing

$$x=\frac{\dot{M}_\mathrm{d}(0)}{\dot{M}_\mathrm{E}}$$
(18)
we obtain

$$\eta \frac{\mathrm{ln}x}{x}$$
(19)
if the mass flow rate decays exponentially, i.e. if (5) is used, and

$$\eta \frac{7}{3}\frac{x^{0.3}1}{xx^{0.3}}$$
(20)
if (4), i.e. (10) holds. In both cases, it is immediately seen that $`\eta `$ is small if $`x1`$, i.e. if during the initial phases of an outburst the mass flow rate is highly super–Eddington.

## 4. The mass flow rate through the outer disc

In order to get $`\eta `$, we need to know $`x`$, i.e. $`\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)`$. Assuming again that $`R_\mathrm{d}=R_\mathrm{h}(\dot{M}_\mathrm{E})<R_{\mathrm{disc}}`$, $`\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)>\dot{M}_\mathrm{E}`$, and using (2), (5), (6) and (8) we find

$$\begin{array}{cc}\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)=\hfill & 3.3\times 10^6\mathrm{M}_{}\mathrm{yr}^1\alpha _\mathrm{h}\alpha _c^{4/5}\hfill \\ & \left(\frac{\mu }{0.6}\right)^1\left(\frac{T_\mathrm{H}}{6500K}\right)(b_1\dot{m}_\mathrm{E})^{1.3}m_1^{0.37}.\hfill \end{array}$$
(21)
We note that $`\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)`$ as given in (21), and therefore also $`\eta `$, still depends on the orbital period. The dependence on $`P`$ enters in a subtle way via the filling factor $``$ which itself must be a function of $`\dot{M}_2`$ and hence via (1) of $`P`$. Unfortunately, this dependence is not explicitly known and can only be determined from full time–dependent calculations of the disc evolution.

Replacing $`\alpha _c^{4/5}`$ in (21) by means of (6), noting that

$$M_d=\dot{M}_2t_\mathrm{q},$$
(22)
we get rid of the poorly known quantities $`\alpha _c`$ and $``$ at the price of introducing $`t_\mathrm{q}`$. This yields

$$\begin{array}{cc}\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)=\hfill & 7.6\times 10^9\mathrm{M}_{}\mathrm{yr}^1(\zeta _e\zeta _R)^1\left(\frac{\mu }{0.6}\right)^1\left(\frac{T_H}{6500K}\right)\hfill \\ & (b_1\dot{m}_\mathrm{E})^{1/4}m_1^{1/2}m_2^{1.7425}p^{0.9281}t_\mathrm{q}(\mathrm{yr}).\hfill \end{array}$$
(23)
Inserting typical values in (23), i.e $`(\zeta _e\zeta _R)1,\alpha _\mathrm{h}0.2,\dot{m}_\mathrm{E}=1.5,m_1=1.4,m_2=1`$ we find $`x1.25t_\mathrm{q}(\mathrm{yr})`$ at $`P=20^\mathrm{d}`$ and $`x10.6t_\mathrm{q}(\mathrm{yr})`$ at $`P=200^\mathrm{d}`$. This shows that if the quiescent time of these discs is long, i.e. $`>`$ a few decades, the accretion efficiency becomes very small. Unfortunately, we do not have any direct observational information about $`t_\mathrm{q}`$. However the analogy with the properties of the outbursts in black hole soft X–ray transients (BHSXTs) (see e.g. King & Ritter 1998) suggests that the quiescent times are likely to be very long, i.e. longer than several decades. In the Discussion below we show that the absence of any observed outbursts from the progenitors of pulsars in wide circular binaries supports this conclusion. We note in passing that $`\dot{M}_\mathrm{d}(R_\mathrm{d})`$ must not be too large if this picture is to be self–consistent: the assumption that $`\nu _{\mathrm{visc}}`$ is dominated by external irradiation, i.e. that we may use (5), requires that

$$\begin{array}{cc}\dot{M}_\mathrm{d}(R_\mathrm{h}(\dot{M}_\mathrm{E}))\hfill & <\frac{8\pi \sigma T_H^4(B_1\dot{M}_\mathrm{E})^{3/2}}{3GM_1}\hfill \\ & <1.2\times 10^5\mathrm{M}_{}\mathrm{yr}^1\left(\frac{T_H}{6500K}\right)^4(b_1\dot{m}_\mathrm{E})^{3/2}m_1^1.\hfill \end{array}$$
(24)
If instead of (5) we use (4), i.e. a viscosity where the disc temperature is self–consistently determined by viscous dissipation, we obtain:

$$\begin{array}{cc}\dot{M}_\mathrm{d}(R_\mathrm{d},t=0)=\hfill & 2.0\times 10^9\mathrm{M}_{}\mathrm{yr}^1(\zeta _e\zeta _R)^{10/7}\alpha _\mathrm{h}^{8/7}\hfill \\ & (b_1\dot{m}_\mathrm{E})^{25/28}m_1^{5/14}m_2^{2.4893}p^{1.3259}t_\mathrm{q}(\mathrm{yr})^{10/7}.\hfill \end{array}$$
(25)
Inserting the same typical values as above in (25) we obtain $`x0.7t_\mathrm{q}(\mathrm{yr})^{10/7}`$ at $`P=20^\mathrm{d}`$ and $`x15t_\mathrm{q}(\mathrm{yr})^{10/7}`$ at $`P=200^\mathrm{d}`$. From this it is seen that if only $`t_\mathrm{q}>`$few yr, disc accretion becomes even more inefficient if the viscosity is due to viscous dissipation rather than external irradiation. We also note that even if (24) is violated initially, the disc’s evolution will eventually be dominated by external heating because $`\dot{M}_\mathrm{d}(R_\mathrm{h})`$ decreases with time.

Finally, in the case of irradiation–dominated viscosity (5) the duty cycle becomes

$$\begin{array}{cc}d\hfill & \frac{\tau _{\mathrm{visc}}(R_\mathrm{h}(\dot{M}_\mathrm{E}))}{t_\mathrm{q}}\mathrm{ln}x\hfill \\ & 0.095\alpha _\mathrm{h}^1\left(\frac{\mu }{0.6}\right)^1\left(\frac{T_H}{6500K}\right)(b_1\dot{m}_\mathrm{E})^{1/4}m_1^{1/2}t_\mathrm{q}^1(\mathrm{yr})\mathrm{ln}x.\hfill \end{array}$$
(26)
Although at a first glance $`d`$ does not seem to depend on $`P`$ it nevertheless does: besides a weak dependence on $`P`$ via $`\mathrm{ln}x`$ there is a significant one from the factor $`t_{\mathrm{q}}^{}{}_{}{}^{1}=\dot{M}_2/M_\mathrm{d}`$. On the one hand $`\dot{M}_2`$ depends explicitly on $`P`$ (cf (1)), on the other hand a more subtle dependence which we have mentioned already earlier enters via the filling factor $``$ in (6). Nevertheless, if typical parameters are used in (26), i.e. $`\alpha _\mathrm{h}0.2`$, $`\dot{m}_\mathrm{E}=1.5`$ and $`m_1=1.4`$, we have $`d0.12t_\mathrm{q}^1(\mathrm{yr})\mathrm{ln}x`$. Therefore, if $`t_\mathrm{q}`$ many decades, $`d1`$.

## 5. Discussion and Conclusions

Li & Wang (1998) have shown that the accretion efficiency of neutron stars which accrete from a giant is very small and, in fact, becomes too small for spinning up the neutron star to millisecond spin periods in systems with a final orbital period $`P_\mathrm{f}>100^\mathrm{d}`$ if the duty cycle $`d<10^2`$. Here we have worked out explicit expressions for $`d`$ and the accretion efficiency $`\eta `$ in terms of the binary parameters and the quiescent time $`t_\mathrm{q}`$ of the outburst cycle in the framework of an irradiated $`\alpha `$–disc model. Our calculations show that both, $`d`$ and $`\eta `$ are small if $`P>P_{\mathrm{crit}}`$ and the mass flow rate at the onset of an outburst is highly super–Eddington. Whether or not $`d`$ and $`\eta `$ are really small enough to prevent neutron stars in systems with $`P_\mathrm{f}>200^\mathrm{d}`$ being spun up to ms–periods depends crucially on the duration of quiescence $`t_\mathrm{q}`$. If $`t_\mathrm{q}>`$ many decades, as some observations of BHSXTs suggest, then both $`\eta `$ and $`d`$ are very small indeed. Unfortunately, determining $`t_\mathrm{q}`$ in the framework of our simple model is not possible. For this fully time–dependent calculations of the disc instability limit cycle would be needed.

However there is another way that we can get a lower limit on $`d`$. We know that pulsars in wide circular binaries must have descended from SXTs. Yet we do not know of a single transient with a period longer than 11.8<sup>d</sup> (GRO J1744–28), even though X–ray satellites would probably have detected an outburst from such a system anywhere in the Galaxy within the last $`30`$ yr. This must mean that outbursts are rather infrequent, i.e. $`t_\mathrm{q}`$ is long. We can make this more precise by defining the following quantities:

Let $`n`$ be the frequency of X–ray outbursts from all systems in the Galaxy, $`N_{\mathrm{PSR},\mathrm{obs}}`$ the currently known number of the pulsar binaries in question, $`N_{\mathrm{PSR}}`$ the current number of the pulsar binaries in question in the Galaxy, $`\tau _{\mathrm{prog}}`$ the lifetime of a progenitor system in the phase of nuclear time scale mass transfer, and $`\tau _{\mathrm{PSR}}`$ the lifetime of the pulsar in the pulsar binaries in question. Here $`\tau _{\mathrm{prog}}`$ and $`\tau _{\mathrm{PSR}}`$ are suitable averages over the relevant population.

Now, the inventory of the binary pulsars in question is incomplete because of flux limitation, dispersion and beaming, whereas an X–ray outburst in one of the progenitor systems would probably be seen wherever in the Galaxy the outburst occurs. We correct for this incompletness by introducing a ‘filling factor’ $`f<1`$. In a stationary situation we then must have

$$t_\mathrm{q}=\frac{N_{\mathrm{PSR}}}{n}\frac{\tau _{\mathrm{prog}}}{\tau _{\mathrm{PSR}}}=\frac{N_{\mathrm{PSR},\mathrm{obs}}}{nf}\frac{\tau _{\mathrm{prog}}}{\tau _{\mathrm{PSR}}}.$$
(27)
$`\tau _{\mathrm{PSR}}`$ must be of the order of the spin-down time scale $`P/2\dot{P}10^8\mathrm{yr}`$, $`\tau _{\mathrm{prog}}`$, on the other hand, is of the order of the nuclear time scale $`t_{\mathrm{}}`$ defined in Ritter (1999, Eq. 14) which, for the systems in question, i.e. those which end their evolution with a long orbital period $`P_{\mathrm{orb}}>200^\mathrm{d}`$, is of the order $`10^8\mathrm{yr}`$. With $`N_{\mathrm{PSR},\mathrm{obs}}=2`$ (i.e. PSR B0820+02 and PSR J1803-2712), $`f0.2`$ and from the fact that we have not seen a single outburst in the past 30 yr this yields the lower limit $`t_\mathrm{q}>(60\mathrm{yr})/f>300\mathrm{yr}`$. We would also like to emphasize that we have adopted a very conservative value for $`f`$. The real value is likely to be smaller still and thus the estimate for $`t_\mathrm{q}`$ even larger because already the canonical beaming factor for pulsars is of order $`0.2`$.

Thus it is very likely that the lack of ms–pulsars in binaries with white dwarf companions and $`P_{\mathrm{orb}}>200^\mathrm{d}`$ is entirely a consequence of dwarf nova–like disc instabilities in long–period binaries.

Finally, we note that spinning up a NS in a long–period binary faces yet another obstacle: during the long quiescent phases the NS acts as a propeller (Illarionov & Sunyaev 1975) and thereby is spun down, making the spin-up to ms spin periods even more difficult.

## 6. Acknowledgment

H.R thanks the Leicester University Astronomy Group, where part of this work was done, for its hospitality and for support from its PPARC Short–Term Visitor grant.

## References

Backer, D.C. Kulkarni, S.R., Heiles, C., Davies, M.M, Goss, W.M. 1982, Nature, 300, 615

Bhattacharya, D., van den Heuvel, E.P.J. 1991, Phys. Rep., 203, 1

Burderi, L., Possenti, A., Colpi, M., Di Salvo, T., D’Amico, N. 1999, ApJ, 519, 285

Frank, J., King, A., Raine, D. 1992, Accretion Power in Astrophysics (Cambridge: Cambridge University Press)

Illarionov, A.F., Sunyaev, R.A. 1975, A&A, 39, 185

King, A.R. 1998, MNRAS, 296, L45

King, A.R., Frank, J., Kolb, U., Ritter, H. 1997, ApJ, 484, 844

King, A.R., Ritter, H. 1998, MNRAS, 293, L42

King, A.R., Ritter, H. 1999, MNRAS, 303, 253

Li, X.-D., Wang, Z.-R. 1998, ApJ, 500, 935

Ludwig, K., Meyer–Hofmeister, E., Ritter, H. 1995, A&A, 290, 473

Ritter, H. 1988, A&A, 202, 93

Ritter, H. 1999, MNRAS, 309, 360

Shakura, N.I., Sunyaev, R.A. 1973, A&A, 24, 337

Taam, R.E., King, A.R., Ritter, H. 2000, ApJ, 541, 329

Tauris, T.M., Savonije, G.J. 1999, A&A, 350, 928

Thorsett, S.E., Chakrabarty, D. 1999, ApJ, 512, 288

van Paradijs, J. 1996, ApJ, 464, L139