Quasielastic scattering in the interaction of ultracold neutrons with a liquid wall and application in a reanalysis of the Mambo I neutron-lifetime experiment

We develop a theory of ultracold and very cold neutron scattering on viscoelastic surface waves up to second-order perturbation theory. The results are applied to reanalyze the 1989 neutron-lifetime experiment using ultracold neutron storage in a Fomblin-coated vessel by Mampe et al. [Phys. Rev. Lett. 63, 593 (1989)]. Inclusion of this theory of the quasielastic scattering process in the data analysis shifts the neutron lifetime value from $887.6\ifmmode\pm\else\textpm\fi{}3$ to $882.5\ifmmode\pm\else\textpm\fi{}2.1$ s.


I. INTRODUCTION
The present reanalysis of our experiment Mambo I [1] of 1989 was stimulated by the current unsatisfactory state of the neutron-lifetime (τ n ) issue.Measurements performed over the past few decades have stated accuracies close to 0.1% [2,3] but disagree by as much as 5 standard deviations.Explicitly, Arzumanov et al. [2] quote τ n = 885.4± 0.9(stat) ± 0.4(sys) s, as compared to the value 878.5 ± 0.7(stat) ± 0.3(sys)s of Serebrov et al. [3].The experiment [1] of Mampe et al. had been the first of a series of fairly high precision measurements using ultracold neutron (UCN) storage in a Fomblin-oil-coated chamber, yielding the value τ n = 887.6 ± 3 s.To separate the small UCN loss due to imperfect wall reflection from that due to β decay, the authors had used a special scaling technique for storage data extrapolation to zero wall loss and applied two corrections to the raw data extrapolation, a gravity correction and a spectral correction.These resulted in a net upward correction of τ n of ∼8 s.More recently, using the theory of quasielastic UCN interaction with viscoelastic surface waves developed since then by Pokotilovski [4] and by Lamoreaux and Golub [5], the data of Ref. [1] have been reanalyzed by Lamoreaux [6] and by Serebrov and Fomin [7].The results were contrary: no correction to the original analysis was proposed in Ref. [6], and a downward shift by ∼6 s for τ n was found in Ref. [7].Our present reanalysis is based on an extended, second-order approach to quasielastic UCN and very cold neutron (VCN) scattering and loss at a liquid wall.This allowed us to include, in the small spectral changes due to slightly inelastic wall collisions, both the transitions from the UCN region to the VCN region, and the reverse transitions from VCNs to UCNs.The latter had been neglected in the previous work [5][6][7].However, both the up and down transitions turn out to be significant for neutrons slightly below or above the critical energy.We take into account a fairly exact model of the trap geometry and experimental conditions of the measurements.The storage lifetimes obtained in our Monte Carlo (MC) simulations with optimized model parameters closely match the values of the experiment, giving us confidence that the model is realistic.
The paper is organized as follows.In Sec.II we describe the basics of neutron quasielastic scattering by thermally induced surface waves of a viscoelastic polymer such as liquid Fomblin oil near room temperature, pointing out the known problems of infinities that occur at very small and/or large values of momentum transfer.Section III introduces our perturbation approach to scattering and the necessity to include second-order perturbation terms to determine the neutron loss probability at a wall reflection.This analysis of reflection, scattering, and loss covers not only the UCN range but also the VCN region where the incident and/or scattered neutrons are in the supercritical range above the total-reflection edge.The full range is required to analyze UCN storage experiments of the Mambo I type where VCNs are not excluded from entering the trap.In Sec.IV we present the analytical results in the form of differential and integral probabilities of scattering back from, and into, the wall, as well as the loss per reflection.We also discuss the numerical methods used to set up the tables used in the Monte Carlo simulations.Section V describes the technical details of the Mambo I experiment, their representation in the computer model, and the strategies used for the simulations.It will be seen that a realistic account of geometry and wall reflection properties is essential to reproduce the 1989 data of measured storage lifetimes versus storage time interval, which are intricately linked to the evolution of the neutron spectrum as a function of storage time.The result of the simulations is the correction needed to determine the neutron lifetime from the measured data.It turns out to be 3.7 s, which is significantly smaller than the correction of 8.8 s originally applied in Ref. [1].The difference appears to be attributable to the slow disappearance of VCNs from the stored neutron spectrum, which occurs since quasielastic upscattering tends to counteract the fast spectral depletion that would be expected from the large wall loss at higher neutron energies alone.This had not been fully taken into account in Ref. [1], although a random walk in energy with small energy changes had been included in the previous analysis.
Applying the new corrections to the uncorrected extrapolation values of the Mambo I experiment we present our revised neutron lifetime evaluation for this experiment in Sec.VI.The discussions in Sec.VII review the current status of the neutron-lifetime value, taking into account also a recent revision of the 2000 double-bottle experiment [2] by the authors of that work.We also discuss the unitarity issue for the Cabibbo-Kobayashi-Maskawa quark mixing matrix, as derived from neutron measurements alone (lifetime and asymmetry parameters).

II. THE DYNAMICAL STRUCTURE FACTOR FOR VISCOELASTIC SURFACE MODES
At high frequencies, the Fomblin oil used in these experiments is expected to behave more like a solid than a liquid, as many polymers do.This behavior can be described by the Maxwell model in the form of a complex-valued, frequency-dependent viscosity where ρ is the density, ν(ω) is the kinematic viscosity, ν 0 is the low-frequency kinematic viscosity, and τ ν is a stress relaxation time characterizing the crossover from viscous behavior at low frequencies to elastic behavior at high frequencies.This model of viscoelasticity was also used in Ref. [4] while purely liquid-like behavior was assumed in Refs.[5][6][7].
Scattering of particles or photons by condensed matter is described in terms of the dynamical structure factor S(q, ω), where q and ω are the momentum and energy transfer, respectively.For thermally excited capillary waves at the surface of an isotropic viscoelastic medium, S(q, ω) has been calculated by Harden et al. [8] in the form where with D 1 (q, ω) = 0 is the surface mode dispersion relation, k B is the Boltzmann constant, T is the temperature, γ is the surface tension (assumed constant), and q and ω are the wave number and frequency of the mode involved in an upscattering or downscattering event.Gravitational effects can be included (as in Ref. [9]) but are neglected here.An alternative form for the capillary wave roughness spectrum S(q, ω) is due to Bouchiat and Meunier [10]: where with Y = γ /4qρν 2 (ω), τ = 1/2ν(ω)q 2 , and s = iωτ .Equation ( 5) had been derived in Ref. [10] for real-valued viscosity (where τ 2 is real-valued and can be taken outside of the argument for the imaginary part), but lengthy algebra shows that the result, written in form (5), is also valid for complex viscosity and that it is equivalent to Eq. ( 2).For calculations we preferred form (2) as it is numerically more stable for large q and ω than form (5).
S(q, ω) is symmetric in q and ω, has a sharp peak at q = 0, ω = 0, and it falls off ∼ω −6 for ω c el q and ∼q −3 for q c −1 el ω, where c el = √ ν 0 /τ ν is the speed of surface wave propagation in the elastic region and ω = c el q is the dispersion relation in this region.
Referring to a normalization issue discussed in Ref. [5], in Eqs. ( 2) and ( 5) we have included factors of 2π in such a way that the usual relation (e.g., Ref. [11]) for static mean-squared roughness amplitude, ξ 2 st = (k B T /2πγ ) ln(q max /q min ), is obtained from S(q, ω): without having to add a factor (2π ) −3 at this stage.Physically, the static surface values, such as ξ 2 st , represent a snapshot of the surface with exposure time t → 0, implying that the range, ω = 1/ t, of the ω integral in Eq. ( 7) is taken over all values from −∞ to +∞.
We note that in quasielastic surface scattering, lower and upper limits on the q range, q min and q max , have to be used to avoid a logarithmic singularity for ξ 2 st and worse than logarithmic for quantities such as the mean-squared surface slope, as shown below.This is due to the power-law behavior of S(q, ω) given in Eqs.(2) or (5).For structure factors with exponential decay, as often used for static roughness, as in Ref. [12], no divergence problem is encountered.It is common practice in experiments on capillary waves, as in Ref. [11], to relate q min to the finite experimental resolution, and q max to the finite intermolecular separation.
The problem of divergences for a fluctuating surface becomes more evident when we attempt to construct a dynamical roughness autocorrelation function of the form where ξ (ρ, t) is the time-dependent surface elevation at point (ρ, 0) = (x, y, 0) of a macroscopically plane surface of large area A, as measured from its average z = 0.The two surface points correlated in Eq. ( 8) are displaced in the plane by δ ρ and in time by δ t .The time integration extends over a long period, t → ∞.
The space-time correlation function C(δ ρ , δ t ) is the Fourier transform of S(q, ω), and the static correlation function C st (δ ρ ) is obtained from C(δ ρ , δ t ) by carrying out the ω integration from −∞ to +∞, for δ t = 0, with the result [5] C Since the integral over the full q range diverges, the series expression in Eq. ( 9) is an expansion, for finite q max and surface points in close spatial proximity, such that δ ρ q max 1.But even with this restriction the leading term, given by Eq. ( 7), remains logarithmically divergent unless restrictions are imposed on both q max and q min .
The dynamical structure factor (2) was used in Ref. [4] in a first analysis of UCN scattering and loss at a viscoelastic wall.In Sec.III we extend this analysis of wall loss to include terms from second-order perturbation theory that are needed for a precise description of prolonged UCN storage as required for a neutron-lifetime measurement.The details are especially important in the critical energy region near the wall potential where the neutrons can cross the UCN/VCN boundary from either side.
The divergence problems alluded to occur in first-order perturbation theory as used in Ref. [4].We might expect an exacerbation if the analysis is extended to second order since, for intermediary states, the integrals over q and ω must be taken over the infinite domain.However, it turns out that no divergences beyond logarithmic occur at this level.
The analysis is very similar to our previous second-order treatment [12] of UCN scattering and loss at a wall with static, rather than dynamic, roughness.The transition from static to dynamic roughness is essentially accomplished by replacing the Fourier transform F (q) of the static height-height correlation in Ref. [12] (see, e.g., the comprehensive review [13], for further details) by the dynamical structure factor S(q, ω) and replacing integrations over q by integrals over q and ω.Since the correspondence in the derivations is almost one-to-one, we will mainly present and apply the final results, referring to Ref. [12] for intermediary steps.

A. Reflection geometry and surface parameters
A plane UCN or VCN wave incident on a planar liquid surface at z = 0 from the vacuum side (z > 0) is of the form (11) with , where m is the neutron mass.We assume that the wave vector k i = (k ix , 0, −k iz ) lies in the (yz) plane.The medium is characterized by the scattering potential where Na is the mean value for the number density N of atoms times their bound-atom scattering length a.The imaginary part η = − Im(q 0 ) Re(q 0 ) describes loss processes (mainly nuclear capture and thermal inelastic scattering).For the per-fluorinated Fomblin oil used in the experiment [1] we can neglect the contribution of incoherent elastic scattering in the supercritical region (k iz > k c0 ), and η is small, of order 10 −5 in the temperature range between 4 • C and 28 • C. k c0 denotes the critical wave number for total reflection at normal incidence on an unperturbed semi-infinite medium with uniform scattering potential.

B. The unperturbed system
For the unperturbed system the complete textbook solution for incident wave (11) with energy hω i , is characterized by the normal wave vector component inside the medium and by amplitudes for reflection and transmission, R i = These results are valid both in the subcritical region k iz < k c0 , where μ i is of order η, and for the supercritical range k iz > k c0 , where μ i is of order 1.At the boundary ), all for η 1.For the incident-wave, flat-wall quantities labeled by subscript i, μ i , R i , and S i , we use the exact expressions, which are valid for any η, to avoid the divergence at k iz = k c0 of low-η approximations such as μ that are commonly used in the UCN region.
Before discussing technical details of the analysis we present a short outline of our approach.
C. Outline of the procedure

Perturbation approach
Using a perturbation approach we expand the coherently reflected and transmitted waves: + ψ (1)  1 (r, t) + ψ (2)  1 (r, t) for z > 0, (13) 1 (r, t) for z < 0, (14) where ψ (1)   1 and ψ (2)  1 , respectively, are the first-order and second-order perturbations propagating into vacuum and ψ (1)   1 and ψ (2)  1 are those propagating into the medium.For energy transfer hω and in-plane momentum transfer hq, the scattered wave vector in vacuum is k whose squared magnitude is The quasielastically scattered intensity is strongly peaked at ω = 0 and q = 0, i.e., in the direction of specular reflection and regular (Snell's law) transmission.[The latter exists as a In-plane momentum transfer q between the projection xk ix of incident wave vector k i and the projection k p of scattered wave vector k.Momentum space is divided into regions where the functions k z and k z in the perturbation integrals are either real valued or pure imaginary.Inside the larger of the two concentric circles, k z > 0; inside the smaller circle, k z > 0; outside the smaller circle, Im k z > 0; and outside the larger circle, Im k z > 0. These boundaries determine the limits of integration over q and ψ.
beam only for supercritical incidence, k iz > k c0 ; for subcritical incidence there is no transmitted beam but the evanescent field carries intensity into the medium with propagation wave number Re(k iz ), which vanishes only for η = 0.] In practice, only scattering events with q and |ω| exceeding certain minimum values q min and , which depend on experimental conditions, are distinguishable from the regular beams.We will return to this issue in Sec.IV.

Approximations
It turns out that the quasielastic intensities are small, of order 10 −4 or less for the Mambo I system.η ∼ 10 −5 is even smaller, and therefore we can use the following approximations for the total normalized intensity reflected and transmitted, for any k ix and any k iz in the UCN or VCN range: for reflection, and for transmission (=loss), where, using symbols as in Ref. [12], I (00) = |R i | 2 is the specular reflectivity, calculated for the given η > 0; I (11) is the intensity per unit of solid angle and energy quasielastically scattered into the vacuum, calculated for η = 0; is the regular transmission, calculated for the given η > 0; I (11) is the intensity per unit of solid angle and energy quasielastically scattered into the medium, calculated for η = 0; I (11) ∼ |ψ (1)  1 | 2 ; I (01) is the interference term between regular reflection and firstorder scattering into vacuum, ψ (1)  1 , calculated for η = 0; I (01) ∼ 2 Re[R i exp(ik iz z) exp(ik ix x) exp(−iω i t)ψ (1) * 1 ]; I (02) is the interference term between regular reflection and secondorder scattering into vacuum, ψ (2)  1 , calculated for η = 0; 1 ]; I (01) is the interference term between regular transmission and firstorder scattering into the wall, ψ (1)  1 , calculated for η = 0; ; and I (02) is the interference term between regular transmission and second-order scattering into the wall, ψ (2)  1 , calculated for η = 0; . The integral intensities scattered into vacuum and wall, respectively, are and p w0 is always less than p D0 and approaches p D0 for k iz k c0 .The two parameters which determine the scattering and wall loss here are the strength of the capillary waves and the absorptive part of the wall potential (η).The surface modes are thermally driven classical waves, as reflected in the dependence S(q, ω) ∼ k B T , in Eqs. ( 2) and ( 5), so we have used k B T as a substitute for ξ 2 , the measure used for static roughness.The above list includes all terms up to order η 1 T 0 and η 0 T 1 in a power-series expansion in η and T , including the constant ∼η 0 T 0 that is nonzero only in the supercritical region.
We have also calculated the terms of order η 1 T 1 but the additional numerical work is significant and the corrections are negligible, as expected.This correction to the loss probability is largest (and negative) at the boundary k iz = k c0 but never exceeds 1% for the oil used in Ref [1].

Loss per reflection
The total loss per reflection, μ(k ix , k iz ), can be calculated in two equivalent ways.The first method, used in Refs.[14,12], is based on subtracting from the incident intensity, which is normalized to unity, the integral reflected intensity (specular plus quasielastic): Alternatively, as in the first-order treatments [4,5], we can add the regular and quasielastic currents crossing into the medium at z = 0 to obtain the integral loss probability If the perturbation approach is taken only to first order, as in Refs.[4] and [5], the last term in Eq. ( 20), I (02) , would be missing.The first-order interference term, I (01) , would still appear but in the subcritical region this term vanishes in the limit η → 0, as shown below.Thus we have agreement with the treatment in Refs.[4,5], which was limited to the subcritical region and based only on the first two terms on the right-hand side of Eq. ( 20).In the supercritical region, I (01) does not vanish, even for η → 0.
The quasielastic loss can occur even for transitions with very small energy changes.Therefore all values of ω down to ω min = −ω i [for real waves with k 2 = 2m h (ω i + ω) > 0] and ω min = −∞ (for virtual states, including k 2 < 0) had to be taken into account.
In the supercritical range we find that the two methods agree perfectly.In the subcritical range k iz < k c0 , method 1 is more stable numerically.The two dominant additions to μ i (k iz ) in Eq. ( 19) are p D0 and I (02) .They are of opposite sign and inspection of Eqs. ( 29) and (33) below shows that both terms are affected, but in the same way, by the vagaries of the choice of q min in the subcritical region where μ i (k iz ) 1, and ≈ 1 and Re k z = 0 in the low-q and low-|ω| range.Adding the two terms in Eq. ( 19) thus removes the uncertainty due to the choice of q min .On the other hand, for method 2 the dominant correction in Eq. ( 20) is p w0 .It also depends on q min logarithmically, and there is no cancellation.We therefore used mainly method 1 and emphasize that the second-order terms I (02) and I (02) are important for either method and cannot be neglected.
After this outline of the procedure, we return to the basics of the perturbation approach, following the method used in Ref. [12], and earlier in Ref. [15].

D. Green's functions
Using the Green's function method, the first-order perturbation ψ (1)  1 in Eqs. ( 13) and ( 14) is obtained as where the integration is performed over the roughness volume V 1 with partly positive and partly negative thickness ξ (ρ, t), with ξ = 0.The integration over t 0 is from −∞ to +∞.This expression is valid for points r on the vacuum side, z > 0, and also for the wave ψ (1)  1 inside the medium, z < 0. The medium is the Fomblin oil film, which is much thicker than the UCN penetration length of ∼20 nm, so it may be approximated as a semi-infinite barrier.The distinction between the vacuum and the medium side appears as a different form of the Green's function G(r, t|r 0 , t 0 ).
As in Ref. [4] and along the lines of [12] we choose a plane-wave representation of the Green's function: The one-dimensional Green's function g(z|z 0 ) for k z takes the following forms [12], with and for z < z 0 , z < 0, and To obtain the second-order perturbation ψ (2)  1 , the zero-order wave function ψ i (r 0 , t 0 ) inside the integral of ( 21) is replaced by ψ (1)  1 (r 0 , t 0 ), and this process could be continued to higherorder perturbations if desired.
Application of this technique to obtain the required quasielastic scattering and loss terms as a function of incident wave numbers k ix and k iz requires straightforward though tedious calculation.We will only present the final results and their application to a reanalysis of the Mambo I experiment.

IV. RESULTS FOR QUASIELASTIC TRANSITION PROBABILITIES
We consider quasielastic scattering from incident UCN or VCN wave vector k i = xk ix − ẑk iz to k = xk x + ŷk y + ẑk z = k p + ẑk z , where we use unit vectors x, ŷ, ẑ with ẑ perpendicular to the wall and pointing away from the liquid surface at z = 0.The in-plane momentum transfer between surface modes and neutrons is hq = h(k p − xk ix ), as shown in Fig. 1, and the positive or negative energy transfer to the Using the technique described in Sec.III, the normalized differential first-order intensity within frequency interval dω, scattered into solid angle d on the vacuum side (z > 0), is obtained in the form with . The last step in Eq. ( 25) follows from the geometrical relationship d 2 q = kk z d since d 2 q is the projection of the area element k 2 d onto the (xy) plane.
The first form for dp D0 in Eq. ( 25) directly corresponds to Eq. ( 20) of Ref. [15] for static roughness scattering (and to similar expressions specifically for the UCN range as in Ref. [12]), except that the structure factor F (q) for static roughness is replaced by the dynamical structure factor S(q, ω), 1  cos θ i ≡ k i k iz is replaced by k k iz , and ω is an additional degree of freedom.
Equation ( 25) can be cast into a form more directly corresponding to Eq. ( 12) of Ref. [4] by using the transition probability w(k i → k) per unit of phase space with where dω = h m k z dk z .Apart from the factors hk z m and 1 2 , Eqs. ( 26) and (12) of Ref. [4] agree.The factor hk z m is, apparently, due to a somewhat different definition of w(k i → k) as used here versus that used in Ref. [4].Regarding the factor 1  2 , we believe that it is needed and that the factor 8 in Eq. (20) of Ref. [4] (which was also used in Refs.[5][6][7]) should be replaced by 4.
It follows from the continuity of the wave function and its gradient at z = 0 that the differential intensity scattered into the wall must satisfy the relation [16] thus the differential transition probability becomes For η = 0, Re(k z ) is nonzero only in the supercritical region The neutron wave propagating into the wall constitutes an enhanced reflection loss in VCN to VCN or UCN to VCN transitions.The next step, required for MC simulations, is the calculation of integral intensities scattered into vacuum (p D0 ) and the wall (p w0 ).It involves triple integration over dω and d 2 q = q dq dψ. Figure 1 shows the projections of incoming and outgoing wave vectors onto the (xy) surface plane and the choice of in-plane angle ψ between − xk ix and momentum transfer q.The geometry is the same as for static roughness (shown in Fig. 1 of Ref. [12]) but, in general, k = k i , and both k i and k can be below or above the critical value k c0 separating the UCN and VCN ranges.Figure 1 shows the case k > k c0 , where k 2 = k 2 − k 2 c0 > 0. For decreasing k the smaller of the concentric circles shrinks and finally disappears at the origin for k < k c0 .
As in Ref. [12], the integrals of ( 25) and (28), and also the integrals needed for I (02) and I (02) , analyzed below, involve terms 2k ix q cos ψ.These integrals over ψ can be performed analytically in terms of elliptic functions, in a way similar to the examples of Appendix D of Ref. [12].The ψ ranges are determined by geometrical requirements (or, for the second-order terms I (02) and I (02) , the absence of limitations in q and ω, as required for intermediary states).For the q value shown in Fig. 1 the ψ range contributing to p D0 extends from −ψ 1 to +ψ 1 , and that contributing to p w0 extends from −ψ 0 to +ψ 0 .The angles ψ 0 and ψ 1 are shown in Fig. 1. k z is real-valued only inside the smaller circle (k p < k ), k z is real-valued inside the larger circle (k p < k), and κ is real-valued everywhere outside the smaller circle (k p > k ).For virtual states, integrations over q and ω are also needed in the range ω < −ω i , where k 2 = (2m/h)(ω i + ω) is negative and both k z and k z are positive imaginary in the entire plane.In this case the ψ integration runs from −π to +π .
We rewrite Eqs. ( 25) and (28) in the form and where with ω max → ∞.Using the analytic expression for the ψ integral, the remaining integrals over q and ω in Eq. (31) were evaluated numerically, usually for a 65 × 65 rectangular grid of points uniform in k ix and nonuniform in k iz , with values from 0 to 4k c0 , for both k ix and k iz .A high density of k iz points was used just below and above k c0 where the total reflection edge is steep though not divergent as long as η > 0. The high density improves interpolation between grid points, as required for the MC simulations.The simulations were performed in almost the same way as in Ref. [12] with the following modification: In calculations of reflection loss from Eqs. (19) or (20), the contributions p D0 and p w0 must include even the smallest positive or negative energy transfers, down to |ω| = 0.The same is true for transitions to and from intermediary states needed to calculate the secondorder interference terms I (0,2) and I (0,2) that will be given below in Eqs.(33) and (35).
However, in determining the next take-off direction and energy in quasielastic scattering back from the wall, in a manner consistent with the probability distribution I (11) (q, ψ, ω), very small energy transfers h|ω| < were excluded from the (q, ψ, ω) search range.For these small energy transfers I (11) is extremely large and the scattered beam is practically indistinguishable from purely elastic reflection or diffraction.More importantly for practical application, the computer search would become forbiddingly slow.As a compromise we usually chose the low-energy limit = 0.3 neV.For the same reason we also had to exclude energy transfers h|ω| > 30 neV from the search.
The search used von Neumann's acceptance and rejection method to determine the scattered-wave vector by use of a set of random numbers; this technique required also a table for the maximum values I max (k iz ) of I (11) in the search range.Except for a very small k iz , the peak position is at hω = at a value q q min that was found numerically.
Applying the methods of Sec.III we obtain for the firstorder and second-order interference terms: For scattering back into vacuum, and where H is the function given in Eq. ( 31).For scattering into the wall, and The first-order terms (32) and (34) contain the static meansquare roughness ξ 2 st given in Eq. (7).Both I (01) and I (01) vanish in the subcritical region where Re(k iz ) → 0 for η → 0. They depend on the choice of q max and q min and cannot be neglected as contributors to the reflection loss in the region k iz slightly exceeding k c0 , from which transitions from VCNs back to UCNs may take place.As a consequence of their signs, both terms enhance the loss.For k iz > k c0 and η = 0, k iz − k iz is real-valued.Thus the last term in the curly braces of Eq. ( 35) vanishes.The integrands in Eq. (31) can be approximated and considerably simplified by setting η = 0.
In the second-order terms, Eqs.(33) and (35), the integral over q in Eq. ( 31) for H (f ) can be extended to infinite range since 2k z falls off sufficiently fast for large |k z | k c0 to ensure convergence.The integral over ω in Eq. ( 31) is also evaluated in the full range, as required for virtual states.We note that an infinite integral of form (31), but with f = 1, appears in the calculation of the static mean-square roughness ξ 2 st given in Eq. (7).For f = 1 the integration over ω can be performed analytically in the complex plane, as in Ref. [5].However, in contrast with the case f = 1, the function f = k z − k z in Eqs.(33) and (35) contains two branch lines in the ω plane.This follows from the ω dependence of k z and k z : . These branch lines add to the branch line for S(q, ω) due to the square root for α(q, ω) in Eq. ( 4).As a result, exact analytical integration in the complex ω plane appeared impossible and we used numerical methods to perform the double integration over q and ω required for I (02) and I (02) .Approximations can be obtained analytically, based on the fact that S(q, ω) is strongly peaked at ω = 0, but those were not used here.
Further simplifications of the interference terms in Eqs.(32)-(35) are obtained, for η → 0, by noting that in the subcritical (supercritical) region, k iz is purely imaginary (real-valued).As a result, I (01) , I (01) , and I (02) vanish for k iz < k c0 .The signs indicate that in the important region just above k c0 , I (02) (I (02) ) makes a negative (positive) contribution to the wall loss probability μ(k ix , k iz ).
We calculated 65 × 65 matrices for the loss probability μ(k ix , k iz ) numerically on the basis of Eq. ( 19), using the same grid as for the integral scattering probabilities p D0 and p w0 .The tables were interpolated in the simulations discussed in Sec.V.
The role quasielastic scattering plays in UCN storage in traps with a liquid wall coating is illustrated in Fig. 2, where we have plotted the mean value μ(k i ) for an isotropic neutron gas.μ(k i ) is obtained from μ(k ix , k iz ) by averaging the incident flux over all angles of incidence.The result is To illustrate the difference between the quasielastic UCN interaction and the elementary flat-wall loss, the plot shows μ(k i ) normalized to the flat-wall isotropic loss probability μ i (k i ).The parameters used for Fig. 2 are listed in the captions; they are characteristic of the conditions of the Mambo I experiment [1]. Figure 2 shows that the loss enhancement is The region where quasielastic scattering from capillary waves is significant extends fairly deep into the subcritical region k i < k c0 .The following parameters were used for these calculations: density ρ = 1897 kg/m 3 , low-frequency kinematic viscosity ν 0 = 2.55 × 10 −4 m2 /s, and surface tension γ = 0.024 N/m as given in Ref. [5], for temperature T = 290 K; viscoelastic transition frequency ν τ = 10 MHz; integration limits q min = 0.8 × 10 −5 nm −1 and q max = 0.8 nm −1 , wall loss coefficient η = 1.65 × 10 −5 , and critical wave number k c0 = 0.07205 nm −1 .pronounced near the critical value k i = k c0 .It extends fairly far into the subcritical region, indicating the importance of quasielastic transitions also outside the immediate vicinity of the critical edge.At the scale of Fig. 2, quasielastic effects are hardly visible above the UCN/VCN threshold but they are not negligible very close to the threshold and could reduce the imbalance in favor of UCN to VCN versus VCN to UCN transitions that is important in UCN storage experiments of the Mambo I type.

V. MONTE CARLO SIMULATIONS FOR MAMBO I
The setup of the Mambo I experiment [1] is shown in Fig. 3. UCNs were confined for variable storage times, up to 3600 s, in a trap of height 30 cm and width 40 cm.The length was variable in the range from ∼1.2 to 55 cm by precise positioning of the vertical rear wall.The vessel was filled and emptied through 65-mm-diameter ports in the front wall, placed at a center height of 8 cm from the bottom.The glass walls were sprayed, and periodically recoated, with liquid Fomblin Y VAC 18/8 oil and kept at a constant temperature in the range from ∼4 • C to 30 • C. The oil served as a low-loss wall material for the UCNs and provided a perfect seal for gaps through which UCNs could otherwise escape.FIG. 3. Setup of Mambo I, the first neutron-lifetime experiment based on UCN storage in a neutron bottle whose walls were coated with liquid Fomblin oil, an excellent UCN reflector.The rear wall of the boxlike vessel is movable to allow variation of the mean free UCN path λ between wall collisions.Most measurements were made with the trap entrance raised 20 cm above the level of the Al window shown on the left, using an S-shaped (rather than straight) stainless steel guide tube connecting the trap to a UCN turbine exit port.Measurement of the storage lifetime as a function of λ and of storage time intervals in Ref. [1] provided the basis for an extrapolation technique for determining the neutron lifetime for β decay.A special scaling technique was used to ensure that the UCN spectrum evolved in the same way for different trap sizes.

A. Features relevant to the simulations
(a) The overall shape deviates from a simple rectangular box, mainly due to sinusoidal undulations of depth 2 mm and wavelength 9-10 mm on the rear wall surface. 1 They were needed to quickly establish isotropy of the UCN gas.As a consequence of the undulations the actual surface area of the rear wall is ∼20% larger than the trap cross section of 30 × 40 cm 2 .This changes the mean free UCN path between wall collisions to a significant extent.
(b) With the sliding shutters closed the entrance and exit openings constitute recessed areas of diameter 65 mm and depth 6 mm.These recesses represent significant deviations from a flat front wall, especially for short trap lengths of 2 cm or less.Both the recesses and the sinusoidal wall profile were taken into account in our simulations.Intersections between the parabolic neutron flight path and the sinusoidal profile were determined in an iterative way, whereas exact solutions were used for the intersections with flat and cylindrical walls.
(c) The UCN storage vessel was connected with a horizontal exit of the UCN turbine source at the Institute Laue-Langevin (ILL) via an S-shaped cylindrical guide with 65-mm inner diameter and a length of 1.1 m.At the turbine exit the UCNs passed through a 0.1-mm-thick Al window that was needed to cut the transmitted UCN spectrum from below.The trap was raised to a level where the entrance port was 20 cm higher, center to center, than the Al window.This geometry ensured that all UCN entering the vessel had enough energy to reach the ceiling in vertical flight, as required for the method of data analysis described in Refs.[1,17].Subsequent to storage the UCNs were guided through a curved section and 65-mm-diameter vertical guide tube to a UCN detector with a 0.1-mm-thick Al entrance window, placed 1.2 m below the exit level.To prevent possible damage to the glass components, both the entrance and exit guide connections left a cylindrical gap ∼1 mm wide.Since these geometrical features are important for the UCN spectrum initially inside the trap after loading, our simulations used albedo, transmission, and loss factors to model the reflection, transmission, and leakage properties of the entrance and exit components and adjust them to the measured emptying time constant of 31 s, for the largest volume.This required an albedo of 0.34 for the exit section leading to the detector.
(d) Closely connected with the issue of reliable geometric parameters for simulations is the question of exact timing.The wall losses depend on UCN energy.Therefore, to ensure virtual independence from the UCN spectrum, which is not known well, and from the spectral changes occurring during storage, scaling Refs.[1,17] was used.For instance, if in Ref. [1] a storage time interval from 1800 to 3600 s was used for a trap length of 55 cm, then for a volume with half the mean free path λ = 4V /S, a storage time interval from 900 to 1800 s had to be used.In this way the UCN spectra developed in the same way for different volumes.The mean free path is evaluated for the actual trap volume V and actual (not projected) surface area S, by taking into account the deviations from a rectangular box [17].All time settings for the end of loading, storage, and emptying were scaled with λ in this way.This implies that the effective storage times, which include the filling and the emptying time constant, are also scaled.
(e) However, both in the experiment [1] and in our simulations, the loading and emptying phases required corrections.In the gas-kinetic approximation the filling and emptying time constants are τ f = τ e = 4V /(S g v ), where S g is the gate cross section.The mean neutron velocity v is almost independent of volume, as shown below (Sec.V E and Fig. 4) by the virtual coincidence of the simulated spectra g1, g2, and g4 for different trap volumes.Thus τ f and τ e scale with volume V , rather than with λ.This estimate neglects the modification (discussed in the next paragraph) due to the filling guide tube section, with a volume of 3.5 l and a decay time of ∼15 s [1], between the Al -foil and the trap entrance, as well as the effect of the exit guide tube between the exit port and the detector.In the experiment [1] and our simulations the mean filling and emptying times were added to the total storage time.
(f) The process of trap loading required further corrections, both for the measurement and for the simulations.Comparing traps with different sizes, we note that there is a small difference in initial spectra, which depends on the characteristics of the filling guide section.The determining factors are the FIG. 4. (Color online) Evolution of neutron energy spectra in the trap as a function of dimensionless energy h h c0 for various storage times.The y axis shows counts per bin for 100 bins over the full energy range.Curve a shows the spectrum, in relative units, when neutrons enter the trap and are launched isotropically, at the entrance port, for path tracking.Curves b to g and i show the spectra at later times.They are the results of MC simulations for the largest bottle volume with a trap length of 55 cm; the spectra for smaller volumes and down-scaled times are practically indiscernible, as shown by the three overlapping curves g1, g2, and g4, for volumes with maximum λ (for g1) and λ reduced by factors of 2 (for g2) and 4 (for g4).Thus the scaling method ensured that the spectra evolve in the same way for different trap volumes.The long-time storage data, for up to 1 h in the experiment and 2 h in the simulations, show that the spectral region around the critical energy h c0 = h2 k 2 c0 /2m (=1.048 m in units of fall height) is depleted very slowly since quasielastic upscattering from the higher-density spectral region below counteracts the large wall loss in this region.This feature plays an important role in determining the neutron lifetime from the experimental data.
following: (i) the velocity-dependent transmission through the Al foil (calculated here as for isotropic incidence on an Al barrier with a uniform scattering potential of 54.1 neV; this neglects the possibility that in the experiment the Al windows may have been contaminated with a thin film of Fomblin oil with higher scattering potential of ∼106 neV due to prolonged pumping of the oil near room temperature although this process is strongly retarded by the low vapor pressure of Fomblin; on the other hand, the difference in mean transmission and effective cutoff energy is reduced for a thin film and by our assumption of isotropic beam incidence instead of the expected forward peaked distribution); (ii) the guide wall loss corresponding to 15-s decay time [1], modeled here as an increase of effective gap width from 1 mm [geometric, as described in (c)] to 2 mm; (iii) the probability of reversal of travel direction for UCNs due to the curvature and surface roughness of the S-shaped feeding guide [described by an albedo factor in (c)]; and (iv) the finite neutron lifetime.β decay during the filling time causes a slight dependence of the initial spectrum in the trap (immediately after filling) since, compared to faster neutrons, the slower ones have a longer filling time constant, and, therefore, are more affected by the β decay.The resulting slight hardening of the spectrum is most significant for the largest volume.In the measurements [1] this deviation apparently was the reason for the fact that the extrapolated endpoint, τ extr , differed by + 1 s ( + 2 s) if derived from the bottle pair with lengths of 55 and 4 cm (1.2 cm), compared to the pair with lengths of 55 and 10 cm.The necessary small correction was applied to the measured storage lifetimes prior to linear regression with error weighting.In our simulations we obtained the same kind of small deviation from linearity, but the deviations were twice as large (1 s → 2 s and 2 s → 4 s).To minimize the uncertainty due to this preregression adjustment we derived our final correction to the neutron lifetime, including the error, using only the simulated data for the three largest volumes.No correction at all was applied to the experimental data for 55 and for 10 cm (and for 55 and 11.2 cm in our simulations), and the straight-line extrapolations, τ extr , for this reference pair served as the common strategy for experiment and simulation.In this way the corrections τ n − τ extr obtained from the simulation for storage time intervals 1800 s/3600 s, 900 s/1800 s, etc., could be directly applied to the experimental extrapolations since these had been obtained on the same basis.(τ n is an estimate of neutron lifetime used as an input parameter for the simulations; see Sec.V B.) In practice, the simulation data extrapolations based on only two points were amended by an error-weighted linear regression that also included the third-largest volume (after the slight spectral correction outlined above).However, the two-point extrapolations practically coincided with those for the three largest volumes.For the neutron lifetime average we obtained from three points τ n = 882.4± 2.1 s, and from two points we obtained τ n = 882.6 ± 2.4 s, for an average result of 882.5 ± 2.1 s.Fits including also the smallest volumes gave consistent results but were not used since, for narrow gaps of only ∼1 cm between the front and the back walls of the trap, they might be skewed by nonuniformity of UCN density over a horizontal plane during filling and emptying.
(g) The trap loading process required yet another adjustment, both for the measurements and the simulations.Since the loading time constant has to be added to the total storage time, the number of wall collisions on the trap walls during filling should also be independent of trap size.To conserve the number of reflections, the aperture of the filling gate was adjusted for each volume.In the simulations a corresponding adjustment of the albedo factor A was made: for instance, from A = 0.18 for 55 cm length to 0.64 (0.67) for trap lengths of 11.2 (4.3) cm, respectively.This corresponds to a reduction of gate area by 44% (40%).
(h) For each pair of holding times, for instance for long holding time t h,long = 3600 s versus short holding time t h,short = 1800 s with the largest trap size, the net storage time, T = t h + τ f + τ e long − t h + τ f + τ e short includes the mean time intervals the neutrons spend in the trap during filling, τ f , and during emptying, τ e .Our simulations showed that the extrapolations were quite sensitive to small variations of T .Underestimating the correction for filling and emptying, τ f + τ e long − τ f + τ e short , by only 1% of τ f + τ e long caused the correction to τ n to increase by ∼1 s.In the experiment the values of τ f and τ e were measured, and in the simulations they are obtained as averages over a large number of trial neutrons (typically 4 × 10 6 runs).
(i) An essential criterion of the reliability of the simulations is a good representation of the measured storage lifetime data τ st as a function of holding times for the pairs 1800 s/3600 s, 900 s/1800 s, etc. (for the largest volume, and scaled down for smaller volumes).To obtain a good fit, we had to determine optimized values for the loss coefficient η (1.65 × 10 −5 for the measurements at 290 K that we simulated) and viscoelastic transition frequency [ν τ = (2πτ ν ) −1 ≈ 10 MHz, where the relaxation time τ ν was defined in Eq. ( 1)].For comparison, τ ν was assumed infinitely large (as for a simple liquid) in Refs.[6,7].Our value of 10 MHz is quite uncertain, especially since the best fit for ν τ depends somewhat on q min and q max , the smallest and largest capillary wave numbers considered.In most simulations we used q min = 0.8 × 10 −5 nm −1 and q max = 0.8 nm −1 .The latter value corresponds to transitions to neutron velocities of ≈50 m/s.We are not aware of any other data for the viscoelastic properties of the Fomblin oil used, but the order of magnitude of megahertz for ν τ is, apparently, typical of other liquid polymers.
(j) Concluding this survey of general considerations for simulations of UCN storage experiments of the Mambo I type, where VCNs were not excluded from the initial spectrum, we emphasize the importance of a reliable model for quasielastic transitions, especially for those from the UCN to VCN range and for return transitions from VCNs to UCNs.In Refs.[6,7] only the former type was considered and a neutron once crossing the UCN to VCN boundary, even if only for the bottom level of the trap (with a gravity boost), was considered lost.In some calculations of Ref. [7] the VCN range was included in the spectra, but only in terms of losses, not quasielastic upscattering or downscattering.In our model both UCN to VCN and VCN to UCN transitions were included.When we ran a "no return" modification, we obtained significant changes in storage lifetimes as a function of storage times.The measured storage lifetime data [1] could not be reproduced, either, when the quasielastic channel was completely turned off, i.e., all reflections were considered purely elastic and specular, or when we used angular scattering distribution, loss, and integral scattering probabilities as for quasielastic scattering but (in violation of the theory) set the energy transfer hω equal to zero.The actual simulations gave a close fit to the experimental storage lifetime data with maximum deviations of <1%.The experimental reference data for the largest volume at 17 • C were τ st = 612.6 ± 3.0 s for storage times 112.5 s/225 s; 639.8 ± 1.3 s for 225 s/450 s; 658.0 ± 1.0 s for 450 s/900 s; 672.0 ± 1.5 s for 900 s/1800 s; and 687.5 ± 2.5 s for 1800 s/3600 s.We do not expect perfect agreement since our model does not include mechanical wall vibrations and imperfect vacuum.In Mambo I wall vibrations were studied by replacing the turbo pump by a vibration-free diffusion pump.The measured lifetimes did not change.The effect of residual gas or vapor on τ n was studied by deteriorating the vacuum from 2 × 10 −7 to 2 × 10 −5 mbar [1].No measurable change of τ n was observed within the statistics of this test.For experiment Mambo II [18], which succeeded Mambo I and used a similar UCN storage system with a liquid Fomblin wall coating, the residual gas composition was measured directly.From these data an upward correction of τ n by <0.1 s is calculated for the stable vacuum of 2 × 10 −7 mbar reported in Ref. [1].We neglect this unidirectional correction since it is insignificant compared to the uncertainty due to vibrations and residual gas that we combine to an estimated net uncertainty of ±0.5 s (Table II).

B. Techniques used for the simulations
Our simulations were performed along the lines described in Ref. [12], with two essential differences: (a) In Ref. [12] all UCNs could be considered to be created at the same start time of a cycle, while in the Mambo I experiment a given UCN could enter the trap at any time between the start and end of the filling period (250 s for the largest volume, both for filling and emptying, i.e., about 8 emptying time constants [1]).Only those entering toward the end of the filling period, within about one loading time constant, had a good chance to remain in the trap until the entrance shutter closed; the rest would be lost or reflected back toward the source.The start time within the filling period for each neutron was determined by a random number.(b) Compared to elastic reflection on walls with static roughness, as in Ref. [12], there is a further degree of freedom, ω.Thus, at each wall impact, a random number determined whether the collision was considered either elastic and specular or quasielastic (with probability p D0 ); if quasielastic, further random numbers determined the values of q, ψ, and ω and thus the initial velocity vector for the next flight segment, in a way consistent with the scattering distribution, which also involves a phase-space factor that slightly favors upscattering over downscattering.
The program keeps track of wall loss accumulation by multiplying a UCN retention factor, initially set equal to 1, by [1 − μ(k ix , k iz )] at each wall collision.At the end of an entire simulated path through the trap the overall loss, including β decay, determines the probability for the neutron to be counted or not.This requires a guess, τ n , as input value for the neutron lifetime.The simulated inverse storage lifetimes, τ −1 st , are obtained from the count rates, plotted versus inverse mean free path λ −1 for the various storage time intervals and extrapolated to λ −1 → 0, as in Refs.[1,17].The result is a value τ extr , and the extrapolation correction τ n = τ n − τ extr represents the sum of gravity correction and spectral correction [1,17].The input value for τ n can be refined and it seems that the correction is quite sensitive to the choice of τ n .We used τ n = 882.4s, which practically coincides with our result τ n = 882.5 ± 2.1 s.However, the change by ∼ + 0.7 s that we obtained for τ n = 881.2s could be an artifact of statistics.The simulations are very slow.

C. Summary of strategy
Our strategy for determining the neutron lifetime τ n from the extrapolation value τ extr may be summarized as follows: The Mambo I system is modeled as closely as possible in terms of all relevant features of the trap itself and its feeding and emptying system, including losses and spectral effects due to the UCN source and modifications by the Al foils used.The model also represents as well as possible the neutron reflection and quasielastic scattering properties of the Fomblin oil used as a wall coating.The simulations yield storage lifetimes τ st as a function of the storage time interval and trap size in the same way as the experiment does.We then choose the same strategy of extrapolating the τ st data to obtain τ extr as was used in the experiment [1] and compare τ extr with the input "guess" τ n (close to the actual neutron lifetime) used for the simulations.The difference τ n − τ extr represents the correction that should be applied to the extrapolation of the τ st data measured in the experiment [1].This method does not rely on whether or not the extrapolation method is "correct" (in terms of directly providing τ n or not) as long as it is the same in the experiment and the simulations.

D. Resulting correction to extrapolations
Our results for the correction are τ n = 5.3 ± 2.1 s for the pair 448 s/896 s of storage times, τ n = 3.7 ± 1.3 s for the pair 896 s/1792 s, and τ n = 2.8 ± 1.0 s for the pair 1792 s/3584 s, where the times apply to the largest volume and are scaled down for smaller volumes.In most simulations we used holding times slightly different from those of the experiment [1], replacing 112.5 s by 112 s, 225 s by 224 s, . . ., and 3600 s by 3586 s.This allowed us to use an integer time step of 1.0 s for all time settings without compromising scaling precision.The difference in τ n between the pairs 1792 s/3586 s etc., versus 1800 s/3600 s etc., is negligible.The final results were obtained with subsecond time steps.They were the same within the statistical error.The corrections for shorter storage times have larger uncertainties and these were omitted.
We used volume settings where λ is 1, 0.5, 0.25, and 0.125 times the value for the largest volume.Even smaller volumes were not used since the narrow gap between the front wall and the rear wall could cause nonuniformity of UCN density over a horizontal surface, as mentioned earlier.On the experimental side, precise measurement of the trap volume and surface area becomes more difficult for small volumes; thus we did not include the data for trap lengths <4 cm in our analysis.For the larger volumes, a 0.1-mm error in measuring the effective trap length (which takes the deviations from a rectangular box into account) would cause a negligible extrapolation error of ∼0.1 s.
In Ref. [1] larger corrections of τ n = 8-9 s were made, based on theoretical analysis [17] for a UCN gas in a boxlike trap under gravity and a spectral correction.The analysis used the elementary, isotropic average of wall loss probability, μ i (k i ), for specular and elastic UCN reflection.Additional simulations also included a small mean energy change, up or down, per collision.It was needed to achieve a good fit to the experimental storage lifetime data as a function of storage times.

E. Role of spectral development during storage
It appears that the details of the quasielastic process at a liquid surface and the resulting spectral development over time, as calculated in this work, are crucial for a more refined description.Figure 4 shows how the UCN/VCN spectra in the trap develop as the storage time progresses, from 112 to 224, 448, 896, 1792, and 3584 s for the largest volume.The spectra include, as curve a, the initial spectrum from the source (inside the trap, i.e., downstream of the Al foil between the turbine source and the S-shaped guide connecting to the trap).Curve a shows a fairly soft Maxwell spectrum, cut from below by the Al window and with an error-function-like upper cutoff of width σ = 0.2 m, centered at energy h = 1.0 m.We use the maximum rise against gravity as a measure of neutron energy, placing the datum at the level of the trap bottom.The soft initial spectrum is consistent with the characteristics of the ILL turbine source [19] test guide that was mainly used for the Mambo I experiment, but the spectrum was not directly measured with this system.The measured storage lifetimes, τ st , were fitted well with this fairly soft initial spectrum.
The sequential storage spectra in Fig. 4 are the simulated spectra, divided by the β-decay factor exp[−(t/τ n )] in order to show more clearly the role of wall loss and spectral change.Curves b to g and i show progressive cooling due to increased wall loss for higher energies, but also by downscattering from regions with higher spectral density.Together with curve g1 for storage time 3586 s we included the corresponding spectra g2 and g4 for mean free paths λ that were smaller by factors of 2 and 4, respectively.The three curves practically coincide.Thus the scaling method works well, ensuring that the spectra evolved in the same way for different trap sizes.
The VCN range is strongly depleted, as expected.However, simulation curve i indicates that even after a storage time of ∼2 h the spectrum is not fully cleaned of neutrons around and slightly above the UCN/VCN boundary h c0 .This region can be replenished by quasielastic upscattering from the higherdensity spectral region below and this tendency is consistent with the simulations of Ref. [7], although these were performed with very different assumptions.Between t ≈ 3600 s and t ≈ 7200 s the spectral density at h = h c0 , referred to its maximum, decreases from 10% to 5%, indicating that cooling eventually prevails over upscattering for this oil.However, for the Mambo I system, the storage times are too long on the scale of the neutron lifetime.Thus, the implicit assumption of analysis [17], namely, the absence of neutrons in the critical region and above, appears violated, and as a result the corrections τ n for the neutron lifetime are smaller than calculated in Refs.[1,17].

VI. RESULT FOR THE NEUTRON LIFETIME
Based on the uncorrected extrapolations of measured storage lifetimes, τ n,uncorr , in Ref. [1] and the new corrections τ n , we obtain the changes to the Mambo I results summarized in Table I.
In Ref. [1] the extrapolation results for the two shortest storage intervals (112.5 s/225 s and 225 s/450 s for the largest volume) had been included in evaluating the overall neutron lifetime.We exclude these data since the uncertainties, both for the uncorrected values from [1] and for the corrections from the present simulations, are much larger and precise scaling becomes difficult, both for the experiment and for simulations.
Averaging the values of Table I we obtain an overall neutron lifetime value of 882.5 ± 1.4 s.Sources of systematic error include minor variations of the correction for different choices of parameters, notably η, q min , q max , τ ν , and the albedo values discussed in Sec.V. We performed test runs for a wide variation of parameters, for instance over two decades each for q min and q max , but the simulations with statistical accuracy sufficient to determine the corrections to the precision presented here were extremely time consuming.The data of Table I and Fig. 4 (for each set of parameters) required about a month of number crunching, employing up to 50 processors.The final results were derived from trap lengths of 55, 11.2, and 4.3 cm, using between 2.4 × 10 7 and 3.6 × 10 7 input neutrons for each trap size.
Table II summarizes systematic effects, both for the experiment [1] and for our MC simulations.The last column contains our best estimates for variations of τ n , up or down, due to the various potential sources of error discussed in this and the previous section.We did not include the slight upward shift of τ n due to residual gas.As discussed in Sec.V A (j), the correction would be <0.1 s for the vacuum conditions reported in Ref. [1] and thus insignificant in view of the net systematic uncertainty of the extrapolations that is estimated at ±1.1 s.The latter value is compatible with the uncertainty of ±1.2 s for the mean value 887.8 s derived from the uncorrected extrapolations for the three longest storage time intervals given in Ref. [1] (last column in Table I).We also excluded from Table II potential uncertainties, both for the experiment and the simulations, due to small detection efficiency variations.If the efficiency for counting following a long storage time is not exactly the same as for the short time, because the spectra are not identical, the storage lifetime τ st derived from this pair is incorrect.However, the effect is very small since for a given trap size τ e changes very little for consecutive storage times, and as long as the experimental system, including the exit system between the trap and the detector, are represented correctly by our model, the deviation will cancel out.The same is true for the fraction of order 10 −3 of UCNs escaping detection, both in the experiment and the simulations, since they are still inside the trap at the end of the emptying period (250 s for the largest volume, and scaled down for smaller trap sizes).
As a result we obtain a systematic uncertainty of ±1.5 s for the correction to τ n .Adding to it the statistical error of 1.4 s in quadrature yields the final result τ n = 882.5 ± 2.1 s.This value and the corresponding values for each of the three time intervals are given in the last row of Table I.We assume that the systematic uncertainty of 1.5 s applies to the final mean value of τ n .
It is difficult to compare these results with the two previous Mambo I simulations [6] and [7], which had also taken quasielastic scattering into account.Both works used the scattering distribution of Ref. [4], which differs by a factor of 2 from our result given in Eqs. ( 25) and (26).Both papers simplified the actual trap geometry with its wall undulations and other details to a rectangular box and both treated quasielastic scattering as a one-way process, where UCNs upscattering to no effect on two-point extrapolations [see Sec.V A(f)] 1 η; q min ; q max ; ν τ ; albedo factors best fit to τ st values strong correlation among η, q min , q max , ν τ 1 Net effect in quadrature 1.5 the VCN range are considered lost.This neglects the possibility of reverse transitions from the VCN range close to the critical energy back to the UCN range.Nor did either paper address the importance of exact timing or the question of whether or not the simulated storage lifetimes as a function of storage time match the data measured in the experiment [1].In our simulations, we included the adjustment of filling gate aperture for each trap size that had been made in the experiment [1].This ensured conservation of the number of wall reflections during trap loading, thus keeping the spectral evolution the same for different trap volumes.Neglecting the adjustment in our simulations resulted in τ n corrections that were ∼1.5 s lower than with the adjustment.While the adjustment was taken into account in [6], no adjustment was reported in Ref. [7].Further differences in the approaches are discussed in Sec.V of the present paper.We believe we have used a more realistic model but refrain from speculating what the overall implications of the simplifications on the correction to an extrapolated neutron lifetime are, or why the authors of [6] obtained a correction consistent with the original correction of 8.8 s in Ref. [1], whereas a smaller average value, 2.8 s, was obtained in Ref. [7].The corresponding mean correction from the present work is 3.7 s.

VII. SUMMARY AND CONCLUSIONS
We have carried calculation of quasielastic UCN and VCN scattering on thermally induced viscoelastic surface waves to second order of perturbation theory and used the results in a new simulation of UCN storage in Mambo I [1], the first dedicated neutron-lifetime experiments based on UCN storage in a liquid Fomblin-oil-coated neutron bottle.Quasielastic upscattering can drive the UCN population over the energy barrier for total reflection, resulting in increased neutron loss from the trap.The second-order corrections are important for a detailed description of spectral development during long storage times of up to 1 h [1], and this behavior critically affects the correction needed to determine the neutron lifetime from the extrapolation of measured storage lifetimes.We obtained smaller corrections and replace the previous result of [1,17] for the neutron lifetime, τ n = 887.6 ± 3 s, by the value τ n = 882.5 ± 2.1 s.The new value is consistent with the 2011 Particle Data Group (PDG) recommendation of τ n = 881.5 ± 1.5 s [20], the average 881.8 ± 1.4 s evaluated in Ref. [18], and the range 880-884 s given in Ref. [21].
In the meantime a new analysis [24] of the UCN doublebottle 2000 experiment [2] resulted in the corrected value τ n = 881.6 ± 0.8 ± 1.9 s.Using the same selection of experiments as in Ref. [20] but replacing the values for Arzumanov et al. [2] and for Mambo I [1] by the new values we obtain the average τ n = 879.8± 0.9 s, where the error is scaled up by the factor of 1.43 derived from χ 2 = 10.2 for five degrees of freedom.The precision stated in Ref. [3] is much higher than for any of the other experiments in the set; as a result, this single measurement carries 63% of the total statistical weight.However, the precision stated in Ref. [3] has been called into question by the authors of Ref. [12], who pointed out a number of possible sources of systematic error for the gravitational trap system that were not taken into account in Ref. [3].Averaging only the five other experiments of the set results in τ n = 882.0± 1.0 s with χ 2 = 2.2 for four degrees of freedom, and thus no error scaling is required.The corresponding confidence level [20] is CL = 0.69, which is significantly higher than the value CL = 0.068 obtained when Ref. [3] is included.
Finally, we examine whether these two neutron-lifetime averages are consistent with the Standard Model (SM) prediction or, possibly, call for an extension at this stage.The up-down quark mixing coefficient |V ud | 2 can be obtained from neutron measurements alone (thus avoiding corrections due to nuclear structure effects) in the form [25] |V ud | 2 = 4908.7±1.9sτ n (1+3λ 2 A ) , where the ratio of axial-vector to vector weak coupling constants, λ A = g A /g V , is determined from polarized neutron decay asymmetry measurements.Using the current PDG average [20], λ A = −1.2701± 0.0025, and τ n = 882.0± 1.0 s, we obtain with the error dominated by λ A uncertainties (which, as in the section on |V ud | in Ref. [20], have been expanded, from 0.0025 to 0.0038, due to experimental inconsistencies).The values (37) and (38) are consistent within 1.0 standard deviations.If we include Ref. [3] in the neutron-lifetime average, and thus use τ n = 879.8± 0.9 s, the divergence between the neutron result, |V ud | 2 = 0.9554 ± 0.0035, and the SM prediction (38) increases to 1.7 standard deviations.
Even with this modest deviation, and more so if we use τ n = 882.0± 1.0 s, it appears that within the current uncertainty level for the neutron measurements (alone), no deviation from the Standard Model needs to be invoked.

.
The reflectivity is |R i | 2 and the loss per unit incident flux is given as

TABLE I .
Results for τ n .

TABLE II .
Summary of systematic uncertainties.