Positron scattering on molecular hydrogen: Analysis of experimental and theoretical uncertainties
Kamil Fedus,1,*Jan Franz,2,†and Grzegorz P. Karwasz1,‡
1Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5/7, 87-100 Toru´n, Poland
2Department of Theoretical Physics and Quantum Informatics, Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, 80-233 Gda´nsk, Poland
(Received 5 March 2015; published 1 June 2015)
Experiments performed in recent years on positron scattering from molecular hydrogen indicated a rise of the total cross section in the limit of zero energy, but essentially disagree on the amplitude of this rise.
Mitroy and collaborators [J.-Y. Zhang et al.,Phys. Rev. Lett. 103,223202(2009)] predicted a scattering length somewhat different from values deduced experimentally. Using a Markov chain Monte Carlo modified effective range theory (MCMC-MERT) we show that the prediction of Mitroy and collaborators allows one to validate the recent experimental results and determine possible uncertainties. By comparing the MERT analysis with the fixed-nuclei density functional calculations we also deduce that probably the effect of virtual positronium formation (or coupling to the virtual positronium state) determines an almost constant value of the total cross section from 3 eV up to the positronium-formation threshold.
DOI:10.1103/PhysRevA.91.062701 PACS number(s): 34.80.Bm
I. INTRODUCTION
The first generation of experiments on positron (e+) scatter- ing from molecular hydrogen (H2) dates back to the 1980s [1–
3]. Those experiments suffered from low counting rates; there- fore they used strong magnetic fields to guide positrons inside the scattering region. Consequently, total cross sections (TCSs) measured in the low energy range were underestimated. In the last decade three new experimental results were published for e++ H2collisions in the low energy range: Karwasz et al.
[4], Zecca et al. [5], and Machacek et al. [6,7]. The latter work also gives differential cross sections (DCSs) measured with the same system. All of the recent experiments report almost twice-higher TCSs in the 3–10 eV range of impact energy when compared to older results. However the agreement between those sets is rather poor, particularly below 1 eV, in spite of the fact that Karwasz et al. [4] and Zecca et al. [5] used essentially the same apparatus [8]. From the theoretical point of view, hydrogen as the simplest representative of molecular targets is a candidate to become a benchmark. Unfortunately, most theoretical results in the low energy range [9–19] are lower than the recent experimental TCSs [4–6].
Mitroy and collaborators using the stochastic variational method, which is one of the most powerful methods for studying few-body systems in the near-zero-energy range [20,21], predicted the mean scattering length of −2.63a0
(here we give a slightly improved value reported recently by Zhang et al. [22]). As they pointed out [20], this result was incompatible with the measurements of Zecca et al.
[5] reported down to 0.1 eV. This conclusion was published before the outcome of the experiment from the Australian National University (ANU) [6] carried out down to 0.5 eV. A good reason to treat Mitroy’s value of the scattering length as correct is the fact that their calculations resolved the
*kamil@fizyka.umk.pl
†jfranz@mif.pg.gda.pl
‡karwasz@fizyka.umk.pl
old discrepancy between theoretical and experimental cross sections for positron annihilation with the H2 molecule (see [23] for more details).
The present work has threefold aims. First, we discuss the advantages and weaknesses of recent experimental systems used to measure TCSs for e++ H2collisions. Secondly, we report a comprehensive study of the compatibility between available experimental TCS and Mitroy and collaborators’
scattering length. For this purpose we carry out a Bayesian predictive analysis of the scattering length using a modified ef- fective range theory (MERT) [24–27] combined with a Markov chain Monte Carlo (MCMC) fitting procedure [28]. Bayesian statistics allows for assessment of the theoretical uncertainties in “classical” MERT analysis, as recently suggested by Editors [29]. In addition we check the compatibility, within the frame of simple phase-shift analysis, between TCSs and DCSs reported by the ANU group [6]. Finally, we report rotationally elastic integral (IECS) and DCSs employing density functional theory (DFT) modeling with the single-center expansion technique (see, e.g., Gianturco and Jain [30]). The calculations are done for energies down to 0.1 eV. Comparing the latter results with MERT analysis we discuss possible discrepancies between theories and experiments occurring in the energy range below the positronium formation threshold.
The paper is organized as follows: In Sec.IIthe comparison between the most recent experimental methods used to study e++ H2collisions is presented. In particular, a critical review of two Trento datasets [4,5] is included. SectionsIIIandIV describe the principles of the DFT and MCMC-MERT models, respectively. SectionVis devoted to the Bayesian predictive analysis of the scattering length. Finally, the paper is concluded in Sec.VI.
II. EXPERIMENTS
In all three recent experiments [4–6] TCSs were evaluated by an absolute method, i.e., measuring the attenuation of the positron beam inside a scattering cell vs the gas pressure.
Karwasz et al. [4] used an apparatus [8] with an electrostatic
[4]
[5]
[6]
[18][17]
[19]
FIG. 1. (Color online) Total cross sections for positron scattering from the ground state of H2. Present DFT results are compared with the calculations of Mukherjee et al. [17], Zhang et al. [18], Zammit et al.[19], and experiments of Karwasz et al. published in 2006 [4]
and the same set presently reanalyzed (Karwasz et al.*), Zecca et al.
[5] and Machacek et al. [6] (directly measured and corrected).
optics and a weak (9–11 G) guiding magnetic field in the region of the scattering cell. The positron source was 22Na salt with 17 mCi activity; the positron moderator was a W-monocrystal 1-µm-thick film, annealed in situ at 3000 K. A compact source- moderator assembly assured a small initial beam divergence (±15◦) and spot (2 mm diameter). A rather long (10 cm) scattering cell with small entrance and exit apertures (1.0 mm diameter) assured a good angular resolution (3.1 × 10−4sr), by more than a factor of 10 better than earlier experiments [1–3].
As discussed in detail for N2[4] using theoretical differential cross sections, a possible underestimation of TCSs at 1 eV in the measurements of Karwasz et al. was as small as 0.4%.
However, the focusing properties of the magnetic field which are highly energy dependent [31] limited the collision energy to the 0–20 eV range.
Data for H2(see normal triangles in Fig.1) published in 2006 [4] were obtained with a rather poor statistics, differently from measurements in Ar and N2, which were repeated with several different settings of optics [4]. Each run consisted of five values of the attenuation ratio I/I0vs a constant pressure, where both I and I0 were mean of 20 values accumulated over 10 s each. The counting rate for I0was between 6 and 9 e+/s in the whole studied energy range (apart from 6 eV where we report the run with I0≈ 4+/s). A limited energy range and few energy points measured in the region 8–10 eV did not allow one to determine precisely the energy bias (i.e., the contact potential + positron work function) between the tungsten moderator and the scattering cell (made of Cu-Ni alloy, ARCAP). For the H2measurements the same bias value (+2.4 eV) was assumed as for Ar and N2 [4], even if we know from the retarding field comparison [32] that this bias could have been by −0.2 eV smaller (but not more than
−0.2 eV). The energy determination is the main source of possible errors in H2measurements since it strongly depends on contact potentials of positron optics (see also Karwasz et al. [33]). A lower pumping speed for hydrogen than for nitrogen while using turbomolecular pumps, makes H2 data particularly vulnerable to the energy bias error and to the
error due to the residual pressure in the scattering cell. To judge on possible systematic errors, H2requires experiments with various settings of the scattering geometry, vacuum arrangements, and beam properties.
Data for H2 from the paper of Karwasz et al. [4] were preliminary results, which did not include whole statistics of measurements performed over the 2004–2005 period in Trento laboratory. In Fig.1(inverted triangles) and in TableI we report TCSs of Karwasz et al. [4] reanalyzed with whole available statistics. TableIincludes also the statistical errors and the number of runs that were performed. The preliminary data published in Ref. [4] are also given for comparison.
Note that data reported in Fig. 7 of the paper by Karwasz et al. [4], includes only few runs at 3–8 eV, giving therefore a somewhat wrong picture of cross sections in this energy range. Additionally the reanalyzed data are reported here assuming +2.2 eV energy bias to show that a small energy shift (within experimental error bars for energy determination) brings Trento experiments to better agreement with the most recent ANU results [6] in the limit of low energy.
Zecca et al. [5] used the same apparatus as Karwasz et al.
[4,8], but by the time of their measurements the radioactive source decayed by almost a factor of 4; the use of a shorter scattering cell (2.1 cm) probably solved some problems of lower beam intensity. However, shortening the scattering cell changes focusing conditions of the magnetic field; in short cells the gas outflow effects need to be tested carefully vs pressure [34]. Differently from Karwasz et al. [4], Zecca et al. applied a +2% theoretical thermal transpiration correction [35] to the measured TCSs; as shown by Blaauw et al. [34], such a correction depends on details of the pumping system and has to be determined experimentally.
The most recent experiment by Buckman and collaborators (ANU Canberra) uses a gas trap [36] to thermalize positrons.
Positrons are emitted from a 50 mCi 22Na source and premod- erated in a solid Ne layer in a backward configuration—a high magnetic field is needed to guide them through the optics and the scattering cell. The Canberra apparatus has two advantages as compared to the Trento setup: a low starting energy of positrons selected for scattering (some 60 meV compared to 1–2 eV from a tungsten moderator) and an intrinsic determination of the collision energy (by the time-of-flight method). However, the use a strong magnetic field (500 G) lowers the angular resolution: at small scattering energies the detector collects all electrons forward scattered. Therefore, a careful correction has to be performed using a retarding-field analyzer (see [37]). Approximate missing angles are up to 33◦ at 0.5 eV and 23◦at 1 eV [38]. The forward-angle correction, say in H2O at 0.5 eV, is estimated as big as 130 ˚A2compared to a TCS of 64 ˚A2directly measured [38]. With similar missing angles the corrections to the H2TCSs are smaller (16% at 1 eV and 6% at 7 eV [6]).
III. THEORETICAL AND COMPUTATIONAL METHODS A. Hamiltonian
In the body-fixed reference frame of the H2molecule, the Schr¨odinger equation is given by
H != Etotal!, (1)
TABLE I. Total cross sections (10−20m2) for positron scattering from molecular hydrogen by Karwasz et al. [4]: preliminary data published in 2006 and the same dataset presently reanalyzed with better statistics and energy shift of −0.2 eV.
Energy (eV) Karwasz et al. [4] Statistical error (no. of runs) Energy (eV) Present reanalysis Statistical error (no. of runs)
1 4.08 0.15 (5) 0.8 4.08 0.15 (5)
1.5 2.705 0.070 (5) 1.3 2.705 0.070 (5)
2 2.144 0.087 (5) 1.8 2.144 0.087 (5)
2.3 1.848 0.068 (5) 2.1 1.848 0.068 (5)
2.5 1.643 0.036 (5) 2.3 1.553 0.074 (11)
3 1.322 0.054 (3) 2.8 1.313 0.050 (5)
3.5 — — 3.3 1.216 0.050 (4)
4 1.219 0.030 (5) 3.8 1.307 0.105 (10)
4.5 1.201 0.032 (5) 4.3 1.216 0.050 (10)
5 1.207 0.015 (3) 4.8 1.255 0.149 (10)
6 1.205 0.019 (2) 5.8 1.097 0.058 (4)
7 1.208 0.018 (2) 6.8 1.063 0.049 (4)
8 1.206 0.022 (2) 7.8 1.118 0.032 (5)
8.5 1.206 0.023 (5) 8.3 1.206 0.023 (5)
9 — — 8.8 1.220 0.050 (5)
10 1.340 0.074 (5) 9.8 1.340 0.0.074 (5)
where Etotalis the total energy and ! is the total wave function.
Within the Born-Oppenheimer approximation and neglecting couplings of the various nuclear and electronic spins, the total Hamiltonian can be written for a fixed nuclear geometry as
H = Te+ Vee+ Ven+ Vnn+ Tp+ Vep+ Vpn, (2) where the various terms are defined below. In all definitions we use atomic units. The kinetic energy of electrons is given by
Te = −1 2
!2 i=1
"i, (3)
where the sum runs over the two electrons. The electron- electron interaction is given by
Vee= + 1
r12, (4)
where r12 is the distance between the two electrons. The electron-nuclei interaction is given by
Ven= −
!2 i=1
!2 A=1
1 riA
, (5)
where riAis the distance between electron i and nucleus A and the sums run over all electrons and nuclei. The nuclear-nuclear interaction is given by
Vnn= +1
R, (6)
where R is the distance between the two nuclei. The kinetic energy of the positron is given by
Tp= −12"p. (7) The interaction between the positron and the molecular electrons is given by
Vep= −
!2 i=1
1 rpi
, (8)
where rpiis the distance between the positron and the electron iand the sum runs over all molecular electrons. The interaction between the positron and the nuclei is given by
Vpn= +
!2 A=1
1 rpA
, (9)
where rpAis the distance between the positron and the nucleus Aand the sum runs over both nuclei.
B. Density functional theory for the positron-molecule interactions
The interaction between the positron and the electrons are described by electron-positron density functional theory. We assume that we can rewrite the Hamiltonian as
H = Hmol+ Hp, (10)
where the molecular Hamiltonian is defined as
Hmol= Te+ Vee+ Ven+ Vnn (11) and Hpdescribes all terms involving the positron
Hp= Tp+ Veff. (12)
The effective potential Veff depends only on the electron density, the nuclear coordinates, and the positron position.
It can be written as
Veff= Vep+ Vpn+ Vpcp. (13) The extra term Vpcp is the electron-positron correlation- polarization potential. The potential Vpcpdescribes the many particle effects between the positron and the electrons and plays a similar role as the exchange-correlation potential in Kohn-Sham theory [39]. Within the DFT description, the total wave function can be written as a product of the electronic wave function #mol(r1,r2) and the positronic wave function ϕ(rp).
!(r1,r2,rp) = #mol(r1,r2) × ϕ(rp). (14)
The electronic wave function is limited to the lowest-energy eigenfunction of the molecular Hamiltonian
Hmol#mol(r1,r2) = EGS#mol(r1,r2) (15) with the eigenvalue EGS. The positronic wave function is the eigenfunction of the effective one-particle Schr¨odinger equation
Hpϕ(rp) = Epϕ(rp), (16) where Ep= Etotal− EGSis the collision energy of the imping- ing positron.
To solve the effective single-particle equation, the single center expansion (SCE) is used for the positron wave function
ϕ(rp) =!
lv
rp−1ulv(rp)Xlv(ˆrp), (17) where Xlvare symmetry-adapted angular basis functions. The indices lv represent the collection of indices (πµlh), where π is the irreducible representation with the component µ and the quantum number l for the angular momentum with component h (for more details see, e.g., Gianturco and Jain [30]). The functions ulv are the corresponding radial parts of the wave function. With this ansatz for the positron wave function, one obtains a set of coupled differential equations for the radial functions
"
d2
drp2−l(l + 1) rp2 + k2
# ulv(rp)
= 2!
l′v′
⟨lv|Veff|l′v′⟩ul′v′(rp), (18) where k =$2Ep and the matrix elements of the effective potential are given by
⟨lv|Veff|l′v′⟩ =
%
d ˆrpX∗lv(ˆrp)Veff(rp)Xl′v′(ˆrp). (19) C. Correlation-polarization potential
The correlation-polarization potential is divided in two parts, depending on the radial distance rp. For large radial distances, the correlation-polarization potential is dominated by the polarization potential
for rp→ ∞ : Vpcp(rp) → Vpol(rp). (20) For short radial distances, the correlation-polarization poten- tial is dominated by the correlation potential
for rp→ 0 : Vpcp(rp) → Vcorr(rp). (21) The two parts are connected at the radial distance rc, which is the outermost point, at which Vpolbecomes larger than Vcorr.
The polarization potential is given by Vpol(rp) = −
"
α0
2rp4 + α2
2rp4P2(cos θp)
#
, (22)
where α0and α2are the values of the isotropic and anisotropic polarizabilities, respectively, and P2(cos θ) is a Legendre polynomial.
The correlation potential Vcorr is modeled by the density functional of Boronski and Nieminen [39]. Boronski and
Nieminen [39] give equations for the energy functional.
Gianturco and Jain [40] give the following explicit equations for the correlation potential in terms of the density parameter:
43rs3ρ−= 1, (23)
where ρ−is the molecular electron density.
For rs! 0.302, Vcorr= 1
2
&
−1.82
√rs + (0.051 ln rs− 0.115) ln rs+ 1.167 '
. (24) For 0.302! rs! 0.56,
Vcorr= 1 2
&
−0.92305 −0.09098 rs2
'
. (25)
For 0.56! rs! 8.0, Vcorr= 1
2
&
− 8.7674rs
(rs+ 2.5)3 −−13.151 + 0.9552rs
(rs+ 2.5)2
− 1.8655
(rs+ 2.5)− 0.6298 '
. (26)
The case rs>8.0 is not considered here, because in this density regime the polarization potential Vpoldominates.
D. Scattering equations
The coupled differential equations in Eq. (18) can be rewritten as coupled Volterra integral equations (see, e.g., Chap. 5 in Gianturco [41] and Franz [42]),
U(R) =% R R0
{J(R) · N(x) + N(R) · J(x)} · V(x) · U(x)dx.
(27) Here J and N are diagonal matrices containing the spherical Ricatti-Bessel and Ricatti-Neumann functions, respectively, V is the potential matrix, and U(R) is a matrix containing radial functions. In the asymptotic region the K matrix can be extracted from the boundary conditions for the wave function.
From the K matrix we can compute the T matrix
T = 1 − (1 − iK) · (1 + iK)−1. (28) Finally the elastic integral cross section can be computed as
σelastic = π k2
!
lv
!
l′v′
((Tllv′v′
((2. (29)
The scattering equations are solved in the body-fixed refer- ence frame. The corresponding K matrices can be transformed into space-fixed matrices by a frame transformation (see, e.g., Chang and Fano [43]). From the space-fixed K matrices we compute rotational elastic and inelastic differential cross sections. In this study, we assume that all molecules are initially in the rotational ground state (J = 0, τ = 0). The final differential cross section is obtained by summation over the initial and final rotational levels.
dσ
d+ =!
J′τ′
dσ
d+(00 → J′τ′). (30)
E. Computational details
The molecular geometry (R = 1.403a0) and electronic structure are computed with the program packageTURBOMOLE
[44]. For all electronic structure calculations the B3LYP (Becke [45], three-parameter, Lee-Yang-Parr [46]) exchange- correlation functional and the aug-cc-pVQZ basis set [47]
are used. The components of the molecular polarizability tensor (α0= 5.545a03and α2= 1.904a30) are computed with time-dependent density functional theory [48]. These values are in good agreement with high-level ab initio values (5.4139a03 and 2.0239a30) of Kolos and Wolniewicz [49] for the lowest rotational-vibrational level. Karplus [50] extracts an experimental value of 5.4367a03for the isotropic polariz- ability from refractive index measurements. Molecular-beam magnetic-resonance measurements by MacAdam and Ramsey [51] of H2 in the J = 1 rotational level gives a value of 2.061 ± 0.003a30for the anisotropy of the static polarizability.
The SCEs of the molecular electron density and of the potential are computed with the computer program package
VOLPOS[52], which is also used to solve the coupled scattering equations by Volterra integration. The grid for the radial integration ranges up to a distance of 146a0. The angular basis set includes terms up to Lmax= 24 for the expansion of the potential and up to lmax= 12 for the expansion of the wave function. At the outermost radial distance the K matrices in the body-fixed frame are generated and elastic integral cross sections σelastic are computed. The transformation of the K matrices from the body-fixed into the space-fixed frame is done with the programPOLYDCS[53], which is also used for the computation of the rotational elastic and inelastic differential cross sections.
IV. MCMC-MERT
A semianalytical approach to MERT, originally introduced by O’Malley et al. [54,55], has been discussed in detail in our previous papers [24,25], so only a brief description will be given here. In particular, we will show how to formulate MERT within the frame of Bayesian statistics.
MERT analytical expression for the partial wave scattering phase shift induced by the spherical part of the long-range polarization potential is given by the following expression (in atomic units) [54]:
tan ηl= m2l − tan δl2+ ˜Bltan δl) m2l − 1* tan δl
)1 − m2l
*+ ˜Bl
)1 − m2ltan2δl
* , (31) where l is the angular momentum quantum number and δl=
π
2(ν − l −12). Here ml and νl denote the energy-dependent parameters which have to be determined numerically using the procedures described in Refs. [24,25]. Parameter ˜Blis related to the additional phase shift that is induced by the unknown short-range potential. This parameter is approximated by the quadratic effective range expansion [24,25]:
B˜l= Bl(0) +12R∗Rlk2+ · · · , (32) where Rl can be interpreted as the effective range for a given partial wave and R∗= α1/2 with α being the dipole polarizability. In the particular case of l = 0, B0 can be
expressed in terms of A, the s-wave scattering length, as B0= −R∗/A.
Using this semianalytical model we verified that in the regime of energies for elastic e++ H2scattering (<10 eV), the leading contributions come from s and p partial waves (l = 0,1) while the contribution of higher partial waves is small and can be described by taking only the leading order contribution to the phase shift [54]:
tan ηl= π αk2 8)
l−12
*)l+12
*)l+32
* for l" 2. (33) Integral elastic (σIE) and differential elastic (dσdω) cross sections are calculated using the single-center partial wave expansions:
σIE= 4π k2
!∞ l=0
(2l + 1) sin2ηl(k), (34)
dσ dω = 1
k2 (( (( (
!∞ l=0
(2l + 1) exp ηlsin ηl(k)Pl(cos θ) (( (( (
2
, (35) where θ is the scattering angle. Substituting Eqs. (31) and (32) for two first partial waves and Eq. (33) for higher angular momentum (up to l = 100 waves in the present analysis) into Eqs. (34) and (35) one gets relations which can be fitted to experimental (or theoretical) cross sections in order to determine the unknown four parameters of the effective range expansions, that is, A(= −R∗/B0), the s-wave scattering length; B1, the zero energy contribution of the p wave; and R0,R1, the effective ranges of two considered waves. All these parameters approximate the contribution of the unknown short-range potential to the scattering phase shift (31).
Unlike in our previous works [24–27] where nonlinear least-square regression procedures were used to determine MERT parameters, here we use a Bayesian inference for parameter estimation [56]. In contrast to the classical fitting, the Bayesian inference does not provide single point estimation in parameter space but rather the probability density functions (PDFs) of model parameters that are shaped by observational (experimental or theoretical) data. Assuming that MERT is correct, we can write the following relation between experimental data and the MERT model:
σexpt.= MERT(µ) + ϵ, (36)
where σexpt.= [σ1,σ2, . . . ,σN] are the experimental data of scattering cross sections (e.g., IECS) and µ = [b0,R0,b1,R1, . . .] is a set (vector) of model parameters appearing in effective range expansion, Eq. (32). Here ϵ = [ϵ1,ϵ2, . . . ,ϵN] denotes a set of random measurement errors plus any real signal in the data that cannot be explained by the model [56]. The starting point for analysis is the Bayes’
rule, which for a multivariate model such as MERT has the following form:
p(µ|σexpt.) = p(σexpt.|µ)p(µ)
p(σexpt.) , (37)
where p(µ|σexpt.) is a joint posterior probability density function of µ after obtaining the experimental data σexpt., while p(µ) is the prior probability density function (which
characterizes knowledge or beliefs about µ before seeing the experimental data σexpt.). Here p(σexpt.|µ) is the likelihood function (the statistical model describing the uncertainties of experimental data) and p(σexpt.) =+
p(σexpt.|µ)p(µ)dµ is a normalization factor, which represents the probability distribution for the experimental data, summed over all possible outcomes µ.
Once a full joint posterior probability distribution is known one can calculate a marginal distribution that is a PDF for any subset of model parameters, i.e., µi ⊂ µ:
p(µi|σexpt.) =
%
p(µi,µ−i|σexpt.)dµ−i, (38) where µ−irepresents the set µ after µiwas removed. Having marginal distributions for each particular parameter it is useful to provide a point estimation representing “best-fit” values together with an estimate of its errors. For example, it can be done using either the mode or the mean value of marginal PDF with the variance of distribution representing its uncertainty [56]. Alternatively, one can give a joint credible region R representing a predictive probability limit of the model due to parameter uncertainties. This can be done calculating the following integral [56]:
%
R
p(µ|σexpt.)dµ = C, (39)
where C is the required credibility (e.g., C = 95%).
Bayesian parameter estimation requires the computation of multidimensional integrals [such as (38) and (39)], and a good solution for this computational problem consists of implementing MCMC methods [56]. Generally speaking MCMC algorithms using prior PDF and likelihood functions generate a sequence {µ(0),µ(1),µ(2), . . .} of model parameters from a Markov chain whose final stationary distribution is a desired p(µ|σexpt.), i.e., the posterior distribution. More details on the background of this sampling numerical method can be found elsewhere [56]. Here we use this computational tool
to estimate MERT parameters in order to describe elastic scattering cross sections for positron collisions with molecular hydrogen. In practice we adapt MCMC MATLAB toolbox by Laine [57] containing the delayed rejection and adaptive Metropolis (DRAM) sampling algorithm with multivariate Gaussian proposal distributions introduced by Haario and co- workers [28]. We assume a Gaussian likelihood p(σexpt.|µ) = N(µ,ν) with the same error variance ν2 for all data points.
The latter is sampled during the MCMC iteration from noninformative conjugate prior inverse 1 distribution (see [57]
for more details). For MERT parameters, we will consider two priors:
(i) Uninformative Gaussian distributions with infinitive variance centered at initial points determined by MERT nonlinear regression fit to experimental data (using the Levenberg-Marquardt nonlinear least-squares algorithm [58]).
Such choice of priors assigns equal probabilities to all possible values of µ components.
(ii) Informative Gaussian distributions for the scattering length with a mean value centered at −2.63a0 and a finite standard deviation of 0.15a0(to include all possible values of Adepending on the geometric position of the molecule with respect to the incoming positron, as reported in [22]).
In all our calculations we use a value of α = 5.314a03[59]
for dipole polarizability. We checked that the final results are practically independent of α in the range of values 5.20a03–5.45a03as reported by different sources (see [59], and references therein).
V. RESULTS AND DISCUSSION
In Fig. 1 we compare presently calculated IECSs using a DFT method with selectively chosen theoretical and ex- perimental data of TCSs and IECSs for e++ H2 collisions below the positronium formation threshold (E < 10 eV).
In the energy range below 2 eV, DFT results are slightly higher than convergent-close-coupling (CCC) calculations of
(a) (b)
FIG. 2. (Color online) MCMC-MERT fits (solid and dashed lines) to directly measured (dots) and corrected (squares) experimental total cross sections by Machacek et al. [6] assuming (a) uninformative priors for all MERT parameters and (b) informative prior for the scattering length with the most probable value of −2.63a0as calculated by Zhang, Mitroy, and co-workers [20–22]. The gray areas in the plot correspond to 50%, 90%, 95%, and 99% posterior regions due to uncertainties of MERT parameters (shown only for the fit to directly measured data). The results of fits are compared with the CCC calculations of Zammit et al. [19] and presently reported DFT.
TABLE II. Posterior mean values of MERT expansion parameters for positron scattering on molecular hydrogen derived on the basis of different experimental data. Here A(= −R∗/B0) is the scattering length, B1is the zero energy contribution for the p wave, and R0and R1 are the s-wave and p-wave effective ranges, respectively. The standard deviations of the final posterior distributions are used as the parameter uncertainties.
A(a0) B1 R0(a0) R1(a0)
Uninformative prior for scattering length
Machacek et al. [6] −2.43 ± 0.03 −2.72 ± 0.30 −1.28 ± 0.10 1.34 ± 0.90
Machacek et al. ang. corrected [6] −2.73 ± 0.03 −2.13 ± 0.07 −0.99 ± 0.03 −1.07 ± 0.28
Karwasz et al. [4] reanalyzed −2.78 ± 0.08 −1.57 ± 0.70 −0.58 ± 0.30 −1.40 ± 1.10
Zecca et al. [5] −2.09 ± 0.02 −4.48 ± 0.90 −1.44 ± 0.10 3.63 ± 1.26
Informative prior for scattering length
Machacek et al. [6] −2.62 ± 0.02 −3.93 ± 1.34 −0.41 ± 0.10 0.15 ± 2.45
Machacek et al. ang. corrected [6] −2.67 ± 0.02 −2.33 ± 0.1 −1.12 ± 0.05 0.18 ± 0.41
Karwasz et al. [4] reanalyzed −2.71 ± 0.02 −0.88 ± 0.10 −0.19 ± 0.20 −1.58 ± 0.44
Zecca et al. [5] No chain convergence achieved
Zammit et al. [19] and rovibrational laboratory frame close- coupling (LFCC) TCS results of Mukherjee et al. [17]. The Rmatrix with pseudostates method by Zhang et al. [18] gives systematically the lowest cross sections in this energy range. At higher energies more discrepancies occur. In particular, DFT becomes lower than CCC data. Moreover, all theories are lower in this energy range (3–10 eV) than the latest experimental TCSs (though the CCC is almost consistent with the reanalyzed data of Karwasz et al. [4] within their experimental error bars).
Forward-angle correction of experimental data cannot explain these differences. The lack of the inclusion of a vibrational inelastic channel in some models also should not be a full explanation since hydrogen is a nonpolar molecule; so direct vibrational and rotational excitation is small, less than 10%
of TCSs according to experimental vibrational cross sections [60] and less than 2% of TCSs according to the LFCC theory [17].
For positron scattering another process is plausible to explain such a difference between theory and experiment, namely, the formation of a virtual positronium, as recently discussed for noble gases by Green et al. [61]. For Ar at 7 eV,
the difference between virtual-positronium [61] and polarized orbital [62] calculations is some 0.7 ˚A2, i.e., of similar value as the difference between the ANU experiment and present DFT for H2at the same energy. Fursa and Bray [63] obtained similar results by modeling coupling to the positronium-formation channel for noble gases. A first step to include this effect within the frame of CCC calculations for e++ H2collisions has been recently done by Zammit et al. [19]. As a result, the latter calculations are higher than other theories (including presently reported DFT) in the energy range of 3–10 eV, as shown in Fig. 1. From a very basic point of view, the short-range interaction such as the virtual-positronium formation can be treated as part of the elastic scattering channel, thus it can be reproduced by MERT parametrization.
Figure 2 shows MCMC-MERT fits to directly measured and corrected TCSs of Machacek et al. [6] using both uninformative and informative Gaussian priors as described in Sec.IV. The correction of the directly measured dataset is related to experimental inability to resolve forward scattering from the incident beam leading to underestimation of TCSs (see Machacek et al. [6] for more details about this correction).
FIG. 3. (Color online) MERT-derived s- and p-wave scattering phase shifts obtained from fits with informative prior to total elastic cross sections by Machacek et al. [6]: (a) directly measured data and (b) corrected data taking into account angular resolution error. For p-wave phase shift only posterior mean values are given. The gray areas in the plot correspond to 50%, 90%, 95%, and 99% posterior regions due to uncertainties of MERT parameters.
[37]
[19] [6]
FIG. 4. (Color online) MERT-derived DCSs obtained from MCMC informative fit to TCSs reported by Machacek et al. [6] (directly measured and corrected datasets). The results are compared with experimental DCS data of Machacek et al. [6] and Sullivan et al. [37] along with theoretical calculations including present DFT and CCC of Zammit et al. [19]. The gray areas in the plot correspond to 50%, 90%, 95%, and 99% posterior regions limits due to uncertainties of MERT parameters (shown only for DCSs derived from directly measured TCSs).
In Fig. 2 MERT-derived IECSs are also compared with present DFT and CCC calculations of Zammit et al. [19].
The uninformative prior provides an excellent agreement between both experimental datasets and the corresponding fits. However, for directly measured TCSs, the fit is lower than both theoretical models (CCC and DFT) below 1 eV, while for corrected TCSs the fit stays in good agreement with presently reported DFT in the low energy limit. On the one hand, the informative prior (assuming Mitory’s scattering length is correct) brings both MERT curves to much better agreement with CCC and DFT at energies below 1 eV.
Moreover for higher energies both curves stay relatively close
to experimental data, nevertheless a fit of much better quality is obtained for the corrected TCS. It indicates that the correction of directly measured TCSs by Machacek et al. [6] is necessary in order to stay in better compatibility with Mitroy’s scattering length.
In Table II we give the mean values and the standard deviations of final posterior distributions obtained from fits to all three recent experimental TCSs [4–6]. Standard deviations are used as uncertainties of theoretical parameters. Except the fit with informative prior to the data of Zecca et al. [5], in general there was no problem with chain convergence (after 16 000 iterations) as indicated by the Geweke test [64] giving
very low |z| values (<1.5) for all MERT parameters. The derived scattering lengths (A) from fits with uninformative prior to presently reanalyzed data of Karwasz et al. [4] and to data of Machacek et al. [6] are slightly higher and lower, respectively, when compared to the mean value of −2.63a0
reported by Zhang and Mitroy [20–22], whereas the value for TCSs of Zecca et al. [5] is noticeably lower. On the other hand, the fits with informative prior to reanalyzed data of Karwasz et al.and Machacek et al. brings the derived scattering lengths to a better agreement with −2.63a0 (still keeping very good agreement between MERT fit and experimental data).
In Fig.3we present the MERT-derived s-wave and p-wave phase shifts from fits to directly measured and corrected TCSs by Machacek et al. [6] using informative priors. For the p-wave only the mean posterior values are presented. The s-wave phase shift derived from directly measured TCSs [Fig.3(a)] predicts the presence of the Ramsauer-Townsend (RT) minimum (a zero value) in the region of 3–6 eV. This is qualitatively similar to the situation in noble gases for Ar, Kr, and Xe [61] where an RT minimum is present at energies preceding the positronium formation threshold (where the cross sections are almost constant). We checked that a similar shape has the s-wave phase shift obtained from presently reanalyzed TCSs of Karwasz et al. [4]; however, predicting the RT minimum with credibility of only 50%. This is due to rather poor statistics of the experimental data [4]. On the other hand the s-wave phase shift derived from corrected ANU TCSs [6] [Fig.3(b)] does not predict the presence of the RT minimum. We note, however, that the correction of directly measured TCSs proposed by Machacek et al. [6] is based on the theoretical differential cross sections by Reid et al. [14] which are three times lower in magnitude than directly measured data. So maybe a more accurate error estimation can slightly change TCSs providing different phase-shift variations than those shown in Fig.3(b). For completeness of analysis we have to add that the energy dependence of s- and p-wave phase shifts derived with uninformative prior (not shown here) have similar shapes as curves in Fig.3(b)with no prediction of RT minimum for all datasets: directly measured and corrected ANU TCSs [6] and reanalyzed data of Karwasz et al. [4].
Finally in Fig.4we check the compatibility between TCS and DCS experimental data reported by the ANU group [6]
(and the San Diego group [37] using a similar technique).
DCSs were calculated using Eq. (35) and the scattering phase shifts shown in Fig.3(we also checked that uninfor- mative prior gives very similar DCSs). Note that in Fig.4we present folded DCSs, i.e., DCS(θ) + DCS(180◦− θ), since experiments are not able to distinguish between forward (θ) and backward (180◦− θ) scattered positrons. For comparison a presently reported DFT calculation and CCC data by Zammit et al. [19] are also shown. Note that both models agree with each other in low energy range, at 0.5 and 1 eV, while for higher energies the discrepancy between them increases.
MERT-derived DCSs from directly measured ANU TCSs are also in very good agreement with the two theories at 0.5 and 1 eV, while DCSs obtained from corrected TCSs are slightly higher. The reason for worse agreement between theories and experiments at these two particular energy points is unknown. Independent measurements are needed in order to confirm experimental DCSs, particularly in such low
energy regime where a noise background is usually high. For higher energies, 3–10 eV, MERT-derived DCSs are in better agreement with experiments than both theories. Moreover, there is a much smaller difference between DCSs derived from directly measured and corrected ANU TCSs despite apparently different phase-shift variations. This in turn validates an almost constant TCS value in-between 3 and 10 eV, as already noticed for several other atomic and molecular targets [65,66]. A hypothesis of the virtual-positronium effects in this energy range is to be exploited theoretically (and experimentally as indicated by the work of Szluinska and Laricchia [67] on noble gases).
VI. CONCLUSIONS
We report calculations of IECSs and DCSs for positron scattering on molecular hydrogen below the positronium formation threshold using a DFT with the single-center expansion technique. The obtained results for IECSs and DCSs are in good agreement in the low energy range (below 1 eV) with recent CCC calculations reported by Zammit et al.[19]. At higher impact energies (3–10 eV), just below the positronium formation threshold, both models diverge.
Moreover, in the same energy range both theories give IECSs lower than experimental TCSs reported recently by Karwasz et al.[4] (presently reanalyzed), Zecca et al. [5], and Machacek et al.[6].
Using MCMC-MERT we carried out a Bayesian pre- dictive analysis of the scattering length using the three experimental TCSs mentioned above. We showed that two datasets, presently reanalyzed Karwasz et al. [4] and Machacek et al. [6], are compatible with the mean value of −2.63a0
obtained with stochastic variational technique by Mitroy and co-workers [20–22]. On the other hand, the TCS by Zecca et al.[5] is incompatible with this scattering length. The same predictive phase-shift analysis shows that DCSs reported by the ANU group [6] are compatible with their experimental TCSs in the energy range between 3 and 10 eV but disagree at lower energies.
The combination of present DFT calculations and MCMC- MERT studies of experimental data indicates that in the low energy range (E < 3 eV) where the experimental TCSs show a rapid rise with decreasing energy, the reanalysis of experiments and/or theories is necessary in order to converge to a final benchmark cross section. The scattering length by Mitroy and co-workers [20] could be a good reference point. In the range of 3–10 eV, where experimental TCSs have almost constant values, a good compatibility between DCSs and TCSs reported by the ANU group [6], as proved by MERT, indicates that the inclusion of virtual-positronium formation inside the elastic scattering channel is probably needed in order to describe total and differential cross sections using ab initio procedures.
Experiments on positron-H2 scattering combined with the tracing of positronium-annihilation processes [67] are needed as well to prove this hypothesis.
ACKNOWLEDGMENT
J.F. thanks Prof. Sigrid Peyerimhoff for making available computational resources.
[1] K. R. Hoffman, M. S. Dababneh, Y.-F. Hsieh, W. E. Kauppila, V.
Pol, J. H. Smart, and T. S. Stein,Phys. Rev. A 25,1393(1982).
[2] M. Charlton, T. C. Griffith, G. R. Heyland, and G. Wright,J.
Phys. B 16,323(1983).
[3] S. Zhou, H. Li, W. E. Kauppila, C. K. Kwan, and T. S. Stein, Phys. Rev. A 55,361(1997).
[4] G. P. Karwasz, D. Pliszka, and R. S. Brusa, Nucl. Instrum.
Methods Phys. Res., Sect. B 247,68(2006).
[5] A. Zecca, L. Chiari, A. Sarkar, K. L. Nixon, and M. J. Brunger, Phys. Rev. A 80,032702(2009).
[6] J. R. Machacek, E. K. Anderson, C. Makochekanwa, S. J.
Buckman, and J. P. Sullivan,Phys. Rev. A 88,042715(2013).
[7] J. R. Machacek, S. J. Buckman, and J. P. Sullivan,Phys. Rev. A 90,042703(2014).
[8] G. Karwasz, R. S. Brusa, M. Barozzi, and A. Zecca,Nuclear Instrum Methods Phys. Res., Sect. B 171,178(2000).
[9] G. Danby and J. Tennyson,J. Phys. B 23,1005(1990).
[10] E. A. G. Armour, D. J. Baker, and M. Plummer,J. Phys. B 23, 3057(1990).
[11] E. A. G. Armour and J. W. Humberston,Phys. Rep. 204,165 (1991).
[12] T. L. Gibson,J. Phys. B 25,1321(1992).
[13] J. L. S. Lino, J. S. E. Germano, E. P. da Silva, and M. A. P. Lima, Phys. Rev. A 58,3502(1998).
[14] D. D. Reid, W. B. Klann, and J. M. Wadehra,Phys. Rev. A 70, 062714(2004).
[15] F. Arretche, R. F. da Costa, S. d. Sanchez, A. N. S. Hisi, E. M. de Oliveira, M. T. d. N. Varella, and M. A. P. Lima,Nucl. Instrum.
Methods Phys. Res., Sect. B 247,13(2006).
[16] S. d’A. Sanchez and M. A. P. Lima,Nucl. Instrum. Methods Phys. Res., Sect. B 266,447(2008).
[17] T. Mukherjee and N. K. Sarkar,J. Phys. B: At., Mol. Opt. Phys.
41,125201(2008).
[18] R. Zhang, K. L. Baluja, J. Franz, and J. Tennyson,J. Phys. B 44, 035203(2011).
[19] M. C. Zammit, D. V. Fursa, and I. Bray,Phys. Rev. A 87, 020701(R)(2013).
[20] J.-Y. Zhang, J. Mitroy, and K. Varga, Phys. Rev. Lett. 103, 223202(2009).
[21] J.-Y. Zhang and J. Mitroy,Phys. Rev. A 83,022711(2011).
[22] J.-Y. Zhang, Y.-J. Yang, Y. Qian, Z.-C. Yan, and U.
Schwingenschl¨ogl,Phys. Rev. A 89,042703(2014).
[23] G. F. Gribakin, J. A. Young, and C. M. Surko,Rev. Mod. Phys.
82,2557(2010).
[24] Z. Idziaszek and G. Karwasz,Phys. Rev. A 73,064701(2006).
[25] K. Fedus, G. P. Karwasz, and Z. Idziaszek,Phys. Rev. A 88, 012704(2013).
[26] K. Fedus,Braz. J. Phys. 44,622(2014).
[27] K. Fedus,Phys. Scr. 89,105401(2014).
[28] H. Haario, M. Laine, A. Mira, and E. Saksman,Stat. Comput.
16,339(2006).
[29] Editors,Phys. Rev. A 83,040001(2011).
[30] F. A. Gianturco and A. Jain,Phys. Rep. 143,347(1986).
[31] A. Zecca, R. S. Brusa, M. Bettonte, E. Rajch, S. Mariazzi, and G. P. Karwasz,Rad. Phys. Chem. J. 68,319(2003).
[32] G. P. Karwasz, R. S. Brusa, and D. Pliszka, Nucl. Instrum.
Methods Phys. Res., Sect. B 251,520(2006).
[33] G. P. Karwasz, R. S. Brusa, and D. Pliszka,J. Phys.: Conf. Ser.
199,012019(2010).
[34] H. J. Blaauw, R. W. Wagenaar, D. H. Barends, and F. J. de Heer, J. Phys. B 13,359(1980).
[35] T. Takaishi and Y. Sensui,Trans. Faraday Soc. 59,2503(1963).
[36] C. Kurz, R. G. Greaves, and C. M. Surko,Phys. Rev. Lett. 77, 2929(1996).
[37] J. P. Sullivan, S. J. Gilbert, J. P. Marler, L. D. Barnes, S. J.
Buckman, and C. M. Surko,Nuc. Instrum. Methods Phys. Res., Sect. B 192,3(2002).
[38] C. Makochekanwa, A. Bankovic, W. Tattersall, A. Jones, P.
Caradonna, D. S. Slaughter, K. Nixon, M. J. Brunger, Z.
Petrovic, J. P. Sullivan, and S. J. Buckman,New J. Phys. 11, 103036(2009).
[39] E. Boronski and R. M. Nieminen,Phys. Rev. B 34,3820(1986).
[40] A. Jain and F. A. Gianturco,J. Phys. B 24,2387(1991).
[41] F. A. Gianturco, The Transfer of Molecular Energies by Col- lision: Recent Quantum Treatments (Springer-Verlag, Berlin, 1979).
[42] J. Franz,J. Math. Phys. 56,012104(2015).
[43] E. S. Chang and U. Fano,Phys. Rev. A 6,173(1972).
[44] Turbomole version 5-6, Copyright ©2002 University of Karl- sruhe.
[45] A. D. Becke,J. Chem. Phys. 98,5648(1993).
[46] C. Lee, W. Yang, and R. G. Parr,Phys. Rev. B 37,785(1988).
[47] T. H. Dunning Jr.,J. Chem. Phys. 90,1007(1989).
[48] F. Furche and R. Ahlrichs,J. Chem. Phys. 117,7433(2002);
,121,12772(2004).
[49] W. Kolos and L. Wolniewicz,J. Chem. Phys. 46,1426(1967).
[50] M. Karplus,J. Chem. Phys. 41,880(1964).
[51] K. B. MacAdam and N. F. Ramsey,Phys. Rev. A 6,898(1972).
[52] R. Curik, F. A. Gianturco, and N. Sanna,J. Phys. B: At., Mol.
Opt. Phys. 33,2705(2000).
[53] N. Sanna and F. A. Gianturco,Comput. Phys. Commun. 114, 142(1998).
[54] T. F. O’Malley, L. Spruch, and L. Rosenberg,J. Math. Phys. 2, 491(1961).
[55] T. F. O’Malley,Phys. Rev. 134,A1188(1964).
[56] P. Gregory, Bayesian Logical Data Analysis for the Physical Sciences(Cambridge University Press, New York, 2005).
[57] Seehttp://helios.fmi.fi/∼lainema/mcmc/.
[58] G. A. F. Seber and C. J. Wild, Nonlinear Regression (Wiley- Interscience, Hoboken, NJ, 2003).
[59] T. N. Olney, N. M. Cann, G. Cooper, and C. E. Brion,Chem.
Phys. 223,59(1997).
[60] J. P. Sullivan, S. J. Gilbert, and C. M. Surko,Phys. Rev. Lett.
86,1494(2001).
[61] D. G. Green, J. A. Ludlow, and G. F. Gribakin,Phys. Rev. A 90, 032712(2014).
[62] R. P. McEachran, A. G. Ryman, and A. D. Stauffer,J. Phys. B 12,1031(1979).
[63] D. V. Fursa and I. Bray,New J. Phys. 14,035002(2012).
[64] S. P. Brooks and G. O. Roberts,Stat. Comput. 8,319(1998).
[65] G. P. Karwasz,Eur. Phys. J. D 35,267(2005).
[66] G. P. Karwasz,Eur. Phys. J. D 37,153(2006).
[67] M. Szluinska and G. Laricchia,Nucl. Instrum. Methods Phys.
Res., Sect. B 221,100(2004).
