Mean volume reflection angle

The mean volume reflection angle of a high-energy charged particle passing through a bent crystal is expressed as an integral involving the effective interplanar potential over a single crystal period. Implications for positively and negatively charged particles, and silicon crystal orientations (110) and (111) are discussed. A generic next-to-leading-order expansion in the ratio $E/R$ of the particle energy $E$ to the crystal bending radius $R$ is given. For positively charged particles, the dependence of the mean volume reflection angle on $E/R$ proves to be approximately linear, whereas for negatively charged particles the linear behaviour is modified by an $E/R$-dependent logarithmic factor. Up-to-date experimental data are confronted with predictions based on commonly used atomic potentials.


I. INTRODUCTION
Volume reflection (VR) is deflection of high-energy charged particles by a planarly oriented bent crystal. It may be regarded as an effect complementary to channeling in bent crystals, in which particles deflect to the side opposite to that of the crystal bending, whereas channeled particles deflect towards the crystal bending. It is considered to be applicable for beam steering at multi-GeV accelerators (e.g., in a multiple-VR mode), having the merit of high acceptance (see [1] and refs. therein). By now, it has been experimentally explored for various bent crystal orientations, incident charged particle types and beam energies [2][3][4][5][6][7][8][9][10][11].
VR effect owes to the asymmetry of the effective continuous potential of atomic planes caused by the crystal bending, but it does not vanish when this asymmetry becomes small. On the contrary, in that case it becomes maximal. Its theoretical treatment is facilitated by the fact that if, as is normally the case, the bending is sufficiently uniform, then the particle angular momentum (or transverse energy) in the effective potential with a centrifugal component is conserved [12]. That makes the problem integrable and in principle analytically tractable [13][14][15][16][17][18][19]. However, since VR angle is accumulated over many interplanar intervals, in each of which the centrifugal potential component is different, yet the interplanar potential itself is generally given by a sophisticated function, the VR angle dependencies on the particle energy and charge sign, crystal material, orientation and bending radius are obscured, ultimately demanding numerical evaluation.
To get a grasp of various parameter dependencies in the VR problem, in [15] it was solved for the simplest example -parabolic model for the interplanar potential. It was demonstrated that the complicated dependencies of the VR angle on R and E greatly simplify if an ex- * Electronic address: bon@kipt.kharkov.ua pansion in the small parameter R c /R [with R c (E) being the critical radius] is carried out. With its aid, explicit expressions for the VR angle of positively and negatively charged particles were obtained. But the model treatment was only adequate for silicon crystal in orientation (110), and did not cover other practically important cases, such as silicon in orientation (111), and germanium crystals.
The aim of the present paper is to generalize the analytic theory of VR, and develop theoretical tools valid for arbitrary interplanar potential and any crystal orientation. To this end, it is expedient to interchange the procedures of radial integration and averaging over transverse energy, which leads to an expression for the VR angle as an integral over a single crystal period of a square root of the effective potential. With its aid, a generic expansion in R c /R can be derived, in which the entire dependence on the straight-crystal potential enters to the coefficient functions. Cases of R ∼ R c and R → R c can be studied, as well.
For simplicity, herein we restrict ourselves to the thickcrystal limit, when contributions from the crystal boundaries are vanishing, and presume the so-called statistical equilibrium, i.e., a uniform distribution in the fast particle transverse energy [13]. Furthermore, we will focus on issues most important for practice, i.e., calculation of the mean VR angle, which is often directly extracted from the experimental data, being almost unaffected by incoherent multiple scattering, and application to silicon crystals, usually used in experiments.
In order to assess the accuracy of the obtained exact and approximate representations, we confront their predictions based on popular realistic potentials with the up-to-date world data for silicon bent crystals in orientations (110) and (111), and for both particle charge signs. For the benefit of the reader, we also include a summary of results for a few model potentials, qualitatively illustrating dependencies on their typical parameters.

II. GENERAL CONSIDERATIONS
When a fast charged particle traverses a bent crystal, the angle θ between its velocity and the family of weakly bent and long atomic planes gradually changes. At its relatively small values, relevant for the VR effect, the particle motion is governed by relativistic classical mechanics in the continuous potential averaged along the planes, similarly to the case of channeling. The period of the force acting on the particle is then tapering to both sides away from the VR region, and in remote regions, in which |θ| ≫ θ c , where θ c = 2V 0 /E is the critical channeling angle for a straight crystal with well depth V 0 , its net deflecting action tends to zero. As for the incoherent multiple scattering, its strength per unit particle path length remains nearly constant everywhere in the crystal (an R-dependent correction to it was evaluated in [16,17]), so, in a thick crystal it eventually becomes formidable. But within the intrinsic VR region, it still remains minor compared to the action of the continuous potential. Besides that, the incoherent multiple scattering is symmetric wrt the initial particle motion direction. Therefore, in the first approximation it should not contribute tot the mean VR angle, and may be neglected.
As was pointed out already in paper [12], where VR was predicted, basic notions facilitating the theoretical description of VR are the same as for channeling in a bent crystal [20]: 1. The continuous potential of a uniformly bent crystal in a planar orientation is axially symmetric. Thus, it conserves the angular momentum projection on the symmetry axis, or, equivalently, the transverse energy including the centrifugal potential. This second integral of the 2-dimensional motion (excluding the irrelevant uniform motion along the crystal bending axis and the planes) makes the problem completely integrable.
2. For a small crystal bending angle (not damaging the ideal crystal lattice), the centrifugal potential may be linearized within the entire crystal volume.
The particle trajectory in polar coordinates {r, ϕ} (radius and angle wrt the bent crystal symmetry axis, located far outside of the crystal) then expresses in an explicit integral form [12]: Here R designates the crystal bending radius, E -the energy of the fast particle (assumed to be ultrarelativistic 1 ), E ⊥ -the particle transverse energy (defined here to include all the r-independent terms), V (r) -the periodic 1 In a generic case, E in Eq. (1) must be replaced by Ev 2 /c 2 , (with period d) continuous potential of the corresponding family of atomic planes, which in a uniformly bent crystal depends only on the radial coordinate r normal to the planes, and r 0 is the r value at the particle entrance to the crystal. Eq. (1) serves as a starting point for all the subsequent calculations. Pure VR experiments, though, measure not the entire particle trajectory inside the crystal, but only the final deflection angle. It can be deduced from (1) as 2 [13] where r c is the radial reflection point, in which + rc R = 0, and factor of 2 in Eq. (2a) accounts for contributions before and after the radial reflection, assuming the crystal to be oriented symmetrically wrt the beam direction. If the first term in the braces in Eq. (2a) is integrated explicitly [15,18], viz., represents the doubled angle between the particle velocity d r/dt (with time differential dt = Rdϕ) and the bent aligned atomic planes at the exit from the crystal, whereas the second term -the polar angle ϕ subtended by the particle trajectory in the crystal, i.e., in effect, the total crystal bending angle. Their difference, naturally, equals the deflection angle in the inertial laboratory frame. Next, it needs to be taken into account that in practice, the initial state represents not a single particle with a perfectly known impact parameter and velocity, but a beam. Particles from the incident beam enter the crystal with random (at an atomic scale) impact parameters and with slightly randomized angles. Accordingly, E ⊥ is where v is the particle velocity and c the speed of light, but insofar as VR is usually applied to ultrarelativistic particles, for brevity we neglect the difference between v and c. For nonultrarelativistic particles, the corresponding substitution is due in the final equations. 2 Compared to [13,15], we define here angle χ of deflection wrt the initial particle direction of motion with the opposite sign, in order to make it positive. Notations θ and θc are reserved for angles of particle motion wrt atomic planes. a randomly distributed variable, too. One has therefore to derive from (2) the angular distribution of the scattering probability. Its calculation is alleviated [13] by relying on the so-called "statistical equilibrium" property, i.e., a uniform distribution in E ⊥ within the relatively small (compared with the initial E ⊥ uncertainty) 3 period ∆E ⊥ = Ed/R: Therewith, the mean VR angle, on which we will focus in the present paper, is given by a simple expression directly in terms of χ(E ⊥ ) rather than its inverse function entering Eq. (3). The constant in the integration limits may be arbitrary, as long as the χ(E ⊥ ) dependence is periodic with the period ∆E ⊥ , and the integration is carried out over the full period.
Although VR is supposed to be formed deeply inside the bent crystal, in principle, integral (2) has yet some residual dependence on its endpoints. The latter, however, fades away with the increase of the target thickness. If, to the leading order in thickness, this dependence is neglected, there remains the genuine volume contribution. Formally, it expresses as a thick-crystal limit: .
is the effective potential of the bent crystal, which has a "washboard" shape due to a tilt introduced by the linearized centrifugal potential component (see Fig. 1). The difference in the braces in Eq. (5) converges on a transverse spatial scale ∆r This longitudinal scale must be smaller than the crystal thickness L = Rθ b , where θ b is the crystal bending angle: The crystals for VR are usually designed to meet the latter requirement. By virtue of the mentioned convergence 3 For a more detailed analysis of conditions of its validity, see [15].
Thick green curve segments mark intervals contributing to the mean VR angle.
at large r 0 , Eq. (5) holds equally well for cases when the VR region is not strictly in the middle of the crystal. For a periodic V (r), the integrand in (5) is not periodic in r. But granted the simplicity of its dependence on the period number, it is feasible to sum over the entire sequence of periods in a closed form [21]. That leads to a generic representation for the VR angle as an integral over a single interplanar period with the integrand given by a generalized Riemann zeta-function [22] (see also [15,17]) with parameter 1/2 and the argument depending on V eff (r): Substitution thereof to Eq. (4) leads to further simplifications, and yields an expression for the mean VR angle involving no special functions: where max is the global maximum in the semiinfinite region r ≤ r ′ < ∞, thus being unique (see Fig.  1), in contrast to possible minor local maxima. Formula (10) can also be derived directly, bypassing evaluation of the angular (or E ⊥ -) distribution, as is demonstrated in Appendix A. The integrand in Eq. (10) differs from zero only in regions V eff (r) < max r ′ ≥r V eff (r ′ ), marked in Fig. 1 by green.
Paradoxically, those are the regions usually associated with channeled particles, whereas VR particles, being over-barrier, on the contrary, radially reflect in the regions marked in Fig. 1 by red. Nonetheless, the latter regions do not contribute to the mean VR angle at all. One of the consequences is that χ must strictly vanish for R ≤ R c (where R c by definition is the smallest value of R, at which V eff (r) has local minima, and hence maxima), insofar as in that case V eff (r) has no maxima. I.e., VR exists under the same condition R > R c (11) as channeling in a bent crystal, but pertains to overbarrier particles. That agrees with the former results of numerical evaluation of integral (2a) in [19] (see also [23]). It is also evident that the integrand of (10) is everywhere positive, so the net deflection always proceeds to the side opposite to that of the crystal bending. As was mentioned in the Introduction, that is the distinguishing feature of VR. In a periodic continuous potential, all of whose wells are equivalent [such as the potential of a silicon crystal in orientation (110)], it is convenient to choose the integration interval end point coinciding with location of any of the maxima r m of the effective potential: max does not depend on r within the entire part of integration interval, where the radicand is greater than zero. Moreover, the integration may be extended over the entire crystal period, provided the region where the radicand turns negative (and the integrand imaginary) is eliminated by taking the real part: That obviates the need for independent maximization of V eff for every value of r in the integrand.
In a more complicated case, when the continuous potential contains two inequivalent wells per period, as it is known to be, e.g., for positively charged particles in a silicon crystal in orientation (111) (see Fig. 5b below), it suffices to split the integration interval in two parts by the intermediate maximum: The subsequent calculation procedure for each of the latter partial integrals is the same as for integral (12) for a single-well potential.
It is also possible [e.g., for negatively charged particles and crystal orientation (111), corresponding to Fig. 5b flipped upside down] that at small E/R the potential is of single-well type, while with the increase of E/R it becomes double-well. In that case, at the critical E/R value one must switch from formula (12) to (13).
Formulas (12), (13) are well suited both for numerical evaluation of the mean VR angle and for analysis of the intrinsic E-and R-dependencies. The latter will be the subject of the next three sections. Phenomenological issues are discussed in Sec. VI.

III. SMALL-Rc/R EXPANSION
For a given oriented crystal, i.e., given V (r), the rhs of Eq. (10) depends on the experimentally changeable parameters E and R only via the ratio E/R, representing the centrifugal force. For any shape of the interplanar potential, the increase of centrifugal force E/R makes the effective potential well progressively more tilted and shallower, wherewith the mean reflection angle (12) decreases. Vice versa, its maximal value is achieved in the formal straight-crystal limit 4 R → ∞. That is the second salient feature of VR mentioned in the Introduction. Since the case of small-R c /R, when χ is maximal, is of the highest practical value, we will analyse it in the first place.
A. Leading order in Rc/R To evaluate the maximal value of χ achieved at R → ∞, one needs merely to replace in Eq. (10) the effective potential by the real one: Since the integrand in the right-hand side of Eq. (14) is a periodic function of r, it does not matter where to choose the integration limits of an interval of the length d.
The obtained result can be cast in a more intuitive form where is the 'local' critical angle (for a straight crystal) sensed by the overbarrier particle during its passage through the interplanar interval with a nearly critical transverse energy, and ... r = 1 d d 0 dr... designates averaging over the continuous potential period. The factor of 2 in Eq. (15) may be interpreted as accounting for contributions from particle motion before and after the reflection point. A corollary from Eq. (15) is that However, in practice, R is never so large compared to R c that the difference between χ and χ 0 is really negligible. Thus, at least a first-order E/R-dependent correction to χ 0 should be taken into account. Its derivation and refinements are discussed below.
B. Next-to-leading order in Rc/R Since Eq. (13) reduces the problem for a two-well periodic continuous potential to that for a single-well one, we begin with the most elementary case when the generic solution is given by formula (12).
To derive from it the next-to-leading order (NLO) correction to (14), note that in Eq. (12) the dependence on E/R enters both to V eff (r) and V eff (r m ), where r m depends on E/R, too. If V (r) has maxima at r = 0 and r = d, the effective potential maximum location r m may be chosen to be close to d. In its vicinity, where we took into account that V ′ (d) = 0. Hence, the potential maximum location shifted by the centrifugal tilt equals This allows to determine the corresponding effective potential maximal value.
which proves to be independent of V ′′ (d) in this approximation. Therefore, to the NLO, i.e., omitting all the O E 2 /R 2 contributions, all the linear E/R dependence in (12) reduces to a term E R (r − d) in the radicand: The change in the integrand compared with the leadingorder (LO) Eq. (14) is illustrated in Fig. 2. Next, linearizing the entire integrand in E/R and integrating the result termwise, we get It involves only the potential for a straight crystal. Since a single-well potential of a straight crystal must be symmetric wrt its minimum at r = d/2, Eq. (20) may be recast as Eq. (21a) demonstrates that the deeper on the average the potential well, or, for a fixed depth, the more rectangular it is, the larger is the LO contribution (the first term), and the smaller by absolute value the negative NLO contribution (the second term). Hence, larger VR angles (beneficial for beam steering applications) are achieved for crystal materials with higher Z and orientations with wider interplanar intervals.
To invariantly relate the χ(E, R) dependence in the NLO with the shape of the potential well, one can introduce a product which, according to the Cauchy-Bunyakovsky-Schwarz inequality, must be greater than unity: for any V (r). The smaller (the closer to unity) this number, the closer the well shape to rectangular.
It should be minded, however, that Eq. (22b) and inequality (23) have been established only for a singlewell potential, and only provided converges, by virtue of which the integration may be extended over the full interplanar interval. Otherwise, notably, for multi-well potentials, the rhs of (22a) may become smaller than unity, because in the limit R c /R → 0 not all the particles completely traverse the last interplanar interval, being reflected in different sub-wells (see examples in Sec. III below).
For a double-well periodic potential, a similar analysis can be based on Eq. (13). Repeating the linearisation procedure for each of the intervals, and presuming the symmetry of each sub-well (which must hold since inversion wrt each sub-well center leaves the periodic sequence of atomic planes invariant), one is led to the result where . .

C. Logarithmic modification of NLO term
Unfortunately, the applicability of simple formula (21a) in practice is undermined because the difference max V − V (r) at the tops of the barriers tends to zero quadratically, wherewith diverges logarithmically. For positively charged particles that happens due to thermal (and zero-point) fluctuations of atom positions in the planes. Granted that those fluctuations are relatively small, a finite result might be attained by just replacing in the NLO integral in Eq. (21a) the interplanar potential by its static counterpart [i.e., letting u = 0 in Eq. (76) below]: with Such an approximation should not be too unreasonable, because if the radial reflection point belongs to the region r c > u, for all r > r c > u the difference between max V stat − V stat (r) and max V − V (r) is small.
Alternatively, one may introduce in (21a) an appropriate logarithmic cutoff r * ≪ d/2: Here, assuming the symmetry of the interplanar potential, we integrated only over half the period and multiplied by 2. The error due to the deliberateness of the choice of r * may be relatively small. For negatively charged particles, however, such an ad hoc cutoff approach may be too crude, because the round top of the potential barrier lies between the planes, being broad for any crystal temperature. In that case, to be rigorous, one should return to representation (19). It can be handled in the next-to-leading logarithmic order by splitting the integration interval in two parts, in the first of which, 0 < r < r i ≪ d, the potential is approximable by its Taylor expansion around the maximum up to the quadratic term: Computing those two integrals and eliminating r i (see Appendix B), we are led to a form similar to (26), with the difference that the lower cutoff is now unambiguously defined and centrifugal-force-dependent: where e = 2.718 is the base of a natural logarithm. The obtained value of r * is natural by the order of magnitude, because the corresponding drop of the potential near the top is commensurable with a typical centrifugal energy Noteworthy, however, is a small numerical factor 1 4e ≈ 0.1. It corroborates our assumption that r * ≪ d/2, as long as Thus, in the case when the particle reflects from the potential barrier in the region of its parabolic top, for description of the NLO correction it suffices to know only two empirical constants: V ′′ (0) and [see Eq. (B2)]. They can be computed numerically, given a realistic parametrization for the interplanar potential V (r). Formula (28) with cutoff (29) is best suited for negatively charged particles, for which the round top of the potential barrier is wide, wherewith the Taylor expansion around its maximum works in a sufficiently broad interval of E/R, too. As for positively charged particles, for them the condition implied at derivation of Eq. (28) may be violated if the tilt proportional to E/R is so large that in the reflection point the continuous potential top is no longer parabolic. But even in that case, around the turnover point the distinction from the extended parabolic approximation for the potential top may be relatively mild (otherwise the NLO approximation itself can break down). Whether or not the logarithmic dependence along with cutoff prescription (29) survive under such conditions will be investigated in more detail on model examples in Sec. IV B 3 and for realistic potentials in Sec. VI.
If the potential has two inequivalent wells per period, with widths d 1 and d 2 , the derivation of the NLO formula must be based on Eq. (13). It gives

IV. MODEL RESULTS
Even though the obtained generic solution (12) or (13) reduces the problem to an integral, which is sufficiently simply calculable numerically, it can be useful sometimes to refer also to model results, in which the dependencies of χ on all the parameters are explicit. Models can also be used for testing the accuracy of the proposed generic approximations, such as the NLO expansion. They may be even more helpful when R becomes commensurable with the critical value R c , so that the NLO approximation breaks down. But it is desirable in this case that the models reflect the shape of the interplanar potential adequately enough. In this section, for references, we will quote predictions for χ for a few such model potentials, pertinent both to single-well and double-well cases.
A. Simple parabolic and square well potentials.
Arbitrary R/Rc Integral (12) can be taken in elementary functions only if V (r) is a linear or quadratic polynomial in r (with an exception considered in Sec. V). In fact, for silicon crystal in orientation (110) the parabolic approximation is known to be rather satisfactory. It is this approximation that was used in [15], but only for small R/R c . More generally, if the effective potential is strongly tilted, only relatively small portions of the wells contribute to the mean VR angle (cf. Fig. 1). Then, if the potential varies smoothly enough, as is typical for moderate-Z crystal materials, it can be parametrized in those small regions by a quadratic polynomial. On the other hand, if the potential walls are steep, as is typical for high-Z materials, one can approximate them by square wells. It will be instructive first to analyse properties of those two simplest among exactly solvable models.
• For a sequence of parabolic wells of depth V 0 > 0, described within a period by the potential serving as an idealization for the interplanar potential for positively charged particles in a silicon crystal in orientation (110), evaluation of integral (12) gives being the critical radius for this model. In the LO, χ 0 = π 2 θ c , while in the NLO, the squared binomial in Eq. (33) linearizes to agreeing with the result of [15], as well as with the generic formula (21a).
Dimensionless parameter (22b) for this model equals • For the sequence of parabolic barriers of height V 0 > 0, being an idealization for the continuous potential for negatively charged particles in silicon in orientation (110), the potential well has the form Inserting this to Eq. (12), one obtains where R c is given by the same Eq. (34). For this model, χ 0 = θ c , while χ0 2d ∂χ ∂(1/R) at R c /R → 0 retains a logarithmic dependence on R c /R. Righthand sides of both (33) and (37) at R → R c tend to zero quadratically, because at R → R c in integral (12) both the width and the height of the integrand Re max V eff − V eff (r) tend to zero linearly. Formula (37) is exact, but follows as well from the generic NLO approximation (28) All the corrections beyond the NLO in this model happen to vanish.
• In the opposite extreme of a square well, the potential within the period has a constant depth V 0 > 0, viz., V (r) = −V 0 within the interval 0 < r < d, and V = 0 at r = 0 and r = d. 5 Then integral (12) evaluates to which for Ed RV0 < 1 matches with the result of [14] in the limit of zero barrier thickness (a → 0). The rhs of this expression as a function of R tends to zero (linearly) only at R → 0, because for a square well R c = 0. The linear law here reflects the fact that at R → R c only the width of the effective potential well shrinks to zero, whereas its depth remains constant, being determined by the height of the sharp walls. Three terms of expansion of (38) read So, for this model, χ 0 = 2θ c [the upper bound of (17) is reached], and [the lower bound of (23) is reached], because θ c (r) here is just a constant. It is also noteworthy that for a square well, the next-to-next-to-leading order (NNLO) contribution to χ is negative, whereas for parabolic well, according to Eq. (33) expanded up to quadratic term in R c /R, it was positive.

B. Double-parabolic potential
More elaborate models may serve to illustrate the dependence of χ 0 and the slope of χ(E/R) on some detail of the shape of the interplanar potential, both for positively and for negatively charged particles. In principle, the integration can be done exactly for any piecewiseparabolic or piecewise-linear potential. For orientation (110), two parabolas can interpolate the planar potential within a period reasonably enough. For orientation (111), where the distances between the planes are not all equal but alternate, producing two inequivalent wells per period (see Fig. 5b below), that is necessary even in the crudest approximation. Let us thus consider the latter case first.

Orientation (111), positively charged particles
The continuous potential for positively charged particles in a silicon crystal in orientation (111) features two unequal potential wells per period (see Fig. 5b below). For simplicity, let us neglect here the thermal smearing of the continuous potential near the planes, treating the potential as consisting of just two alternating parabolic wells of different widths d S , d L and depths V S , V L > 0. To facilitate correspondence with the cases of negatively charged particles, as well as orientation (110) considered below, it is expedient to write the continuous potential within a period in a symmetric form: where and are the halves of the two deep and broad wells, and with d L < d, V S < V L , describes a minor midway well. Under such conditions, we can apply formula (13). Evaluation of the corresponding partial integrals gives At that, Therefore, the mean VR angle (45) in this case amounts the average of critical angles for each of the sub-wells [cf. Eq. (33)], with weights proportional to the widths of the corresponding intervals. For the limiting value χ 0 that is rather obvious, in view of representation (15). Its validity through all orders in R c /R owes to the similarity between the two unequal wells.
It illustrates that inequality (23) may break down for a double well.

Orientation (111), negatively charged particles
For negatively charged particles in the same crystal orientation (111), the continuous potential is inverted upside down. Then, in effect, it contains just one well, which merely features a small bump in the middle of its bottom. For the present case, calculation by exact formula (12) leads to a rather bulky result, so, for simplicity, we will restrict ourselves to the NLO calculation.
The corresponding limiting angle (14) is computed as If V S < V L /2, arth in Eq. (50) may be expanded as Putting there values d S /d = 1/4 and V S /V L ≈ 1/3 corresponding to a real silicon crystal, one finds That is somewhat lower than for positively charged particles [cf. Eq. (48)], although the difference is not as significant as in the case of orientation (110) [cf. Eqs. (33) and (37) at R/R c → ∞].
The slope of the √ Eχ dependence on E/R can readily be evaluated, too. In total, Eq. (53) is valid for any values of V L , V S , d L , d S , provided the coefficient at arth is positive: Physically, that implies that the particle must have enough transverse energy to overpass the midway bump in the last interplanar interval of the bent crystal, in spite of the centrifugal tilt. It is also straightforward to check that in the formal limit d S → 0, the second line in (53) vanishes, and the result boils down to Eq. (37) for a sequence of purely parabolic barriers [relevant for orientation (110)]. In the opposite formal limit d L → 0 (a narrow major barrier), and V S → 0 (a wide but shallow bump), (53) goes over to the result (39) for the square well to the NLO in R c /R [which was written there for positively charged particles in a high-Z crystal in orientation (110)].

Orientation (110), positively charged particles
The pure parabolic model for orientation (110) considered in Sec. IV A does not take into account the smoothness of the continuous potential in vicinities of atomic planes (cf. Fig. 5a below), which, as was shown in Sec. II, can give rise to a logarithmic modification of the linear E/R dependence. This shortcoming can be amended by adopting a double parabolic potential, among the pieces of which, in contrast to the situation of Secs. IV B 1, IV B 2, one is convex upwards, while the other is convex downwards: Here is the harmonic bottom of the potential well, with d 2 < d and V 2 < V 0 , whereas parts The adjustable parameters in this model are the total well depth V 0 and the shape parameter Potential (55) may be substituted to Eq. (12). As long as V 1 is small, in practice it will normally satisfy the inequality wherewith the reflection point will belong to the domain of V 2 rather than V 1 [see Eq. (55)]. Condition (61) can be solved for R as where the right-hand side is much greater than provided d 1 ≪ d 2 . Since condition (62) is usually met in practice, strictly speaking, in this model R can not be regarded as asymptotically large. Therefore, the NLO formula does not strictly apply here, and to be rigorous, one has to return to Eq. (12). In Appendix C it is shown that a satisfactory approximation for χ under the present conditions is Compared to Eq. (33), the main correction here stems from the overall factor It also appears that in spite of the smearing of the potential around the planes, the dependence of (64) on E/R under condition (62) does not contain a logarithmic factor, because the particle merely does not enter the smeared region at all. The model (55) is already sufficiently detailed to be confronted with experimental data [5,9,24] available for silicon crystal in this orientation. Its fit used in the NLO equation (28) gives parameter values It is also close, with nearly the same but less tightly constrained parameters, to the calculation by exact formula (12) for model (55).

Orientation (110), negatively charged particles
To obtain from (55) the potential for negatively charged particles, there is no need to flip it and shift by half the period -it suffices just to replace δ by 1−δ, which interchanges the small and large parameters: d 1 ↔ d 2 , V 1 ↔ V 2 . Since a negatively charged particle traverses all the three regions of piecewise-parabolic potential (55), the result of the exact calculation is more cumbersome. But the NLO calculation is simple. Evaluation by Eqs. (14) and (55)-(60) gives with θ c = 2(V 1 + V 2 )/E and δ = d 1 /d ≪ 1, while evaluation of the slope can be done straightforwardly by Eq. (28): That explicates the δ-dependent correction to Eq. (37). Naturally, the increase of δ in turn increases χ 0 [in contrast to the case of positively charged particles -cf. Eq. (64)], whereas the slope diminishes, because the potential well broadens [see the commentary after Eq. (21)].
The expressions for χ obtained in this subsection for double-parabolic model potentials exhibit surprisingly simple dependencies on various typical parameters. They can even provide rather accurate phenomenological predictions, which can be used for guidance in practice. But the meaning of such coincidences should not be overestimated.

V. R → Rc LIMIT
It may also be of interest to determine the behaviour of χ near the critical point R = R c , although it can hardly be of high practical value, insofar as χ there tends to zero. In the previous section we had seen that for parabolic potential models, χ vanishes at the threshold quadratically: But those models presume that the second derivative of the effective potential remains constant over the entire interplanar interval. In reality, the potential is smooth near the atomic planes. At R → R c the potential well is tilted so strongly that the radicand of (12) is positive only in a small vicinity of point r s , wherein V ′ (r s ) is maximal, i.e., V ′′ (r s ) = 0. To describe the effective potential under such a condition, one thus has to expand V (r) up to the third order in r − r s : Here V ′ (r s ) = F max > 0, and V ′′′ (r s ) < 0. Therewith, integral (12) assumes the form 6 This relation is valid both for positively and for negatively charged particles. Compared to (69), the decrease law in (71) is steeper, but the prefactor involves |V ′′′ (r s )| −3/4 ∼ u 3/2 , with u being the rms thermal displacement of atoms off the plane [see Eqs. (76) and (77) below], whereby ratio is suppressed by a small factor (u/d) 3/2 . Therefore, asymptotic law (71) holds only in a narrow vicinity of R c , whereas at R ∼ 4R c the behaviour of χ/θ c may be more similar to (69).

VI. USE OF REALISTIC POTENTIALS
Let us finally turn to a more realistic description of the mean VR angle, based on exact generic representations (12), (13) for χ. Our knowledge of the interplanar potential, even though good, is in principle only approximate, and its uncertainty gives rise to some uncertainties in theoretical predictions. We will complete our study by assessing this sensitivity. That can be done best in comparison of predictions obtained with different potentials with experimental data. A comparison of experimental data with the potential deflection theory was formerly presented in [19] for positively charged particles. Later, there were also published detailed data for negatively charged particles reflecting on a crystal in orientation (111) [10]. Along with the phenomenological check of exact equations, it is important yet to assess the accuracy of NLO approximations derived in Sec. III.

A. World data. Intercepts and slopes
Experimental data on VR, obtained for various crystals and their orientations, for various incident particles and their energies, can be combined together by plotting √ Eχ vs. E/R. Up-to-date data (see Figs. 6-8 below) do not exhibit a marked curvature, basically being compatible with linear dependencies. Thus, for those cases the logarithmic modification of the slope calculated in Sec. III and present in theoretical predictions for negatively charged particles in Figs. 8, 9 is not revealed yet, because of the scatter between the points. Nonetheless, two numbers inferred from experimental data (the intercept and the slope of the linear dependence) for each particle charge sign and crystal orientation together can appreciably constrain the interplanar potential.
The tightest experimental constraint is for χ 0 , it corresponds to an accuracy of a few percent (standard deviation) -see Table I below. The slope d(E/R) is constrained more loosely -only to within 5-30%. Nonetheless, it can already rule out the crudest models, such as the pure parabolic potential for (110) orientation, containing only one parameter. For instance, for 400 GeV protons on Si (110) [5,9,24] from the linear regression to the experimental data one infers which is more than 6 standard deviations below the corresponding pure parabolic potential model value (36). It must be realized, though, that the experimentally covered range of curvatures extends somewhat beyond 1/4R c , wherewith the NLO approximation itself becomes insufficiently accurate, and the value of K determined in such a broad interval of 1/R appears to be underestimated due to sizable NNLO corrections. The use of Eq. (33) valid beyond the NLO improves the agreement, but does not make it perfect, anyway. A double parabolic potential model of Sec. IV B 3 has two adjustable parameters, by virtue of which it can sustain statistical comparison with VR experiments better.

B. Parametrizations for atomic potentials
Most valuable, of course, would be to compare the experimental data with predictions based on the best available potential. There are ab initio calculations thereof at a varying degree of sophistication, but the crystal structure is in principle more complicated than that of its constituent atoms. Therefore, a commonly adopted attitude for construction of crystal potentials is just to add potentials for isolated atoms, with Z 1 being the projectile charge in units of proton charge e, by placing them at the lattice sites, and neglecting interaction between the atoms, i.e., solid-state binding effects. When the number of covalently paired electrons is much smaller than the total electron number, their effect is expected to be relatively small. For silicon, there are 4 valence vs. total 14 electrons (although half of the covalent bonds is extended along the considered low-index planes), so, the validity of this approximation may be not self-evident . Nonetheless, in this paper we confine our objective to comparing predictions from noninteracting-atom potentials. The potential of a thus obtained atomic plane has the form where n s is the areal density of atoms in the plane. It is also necessary to average it over the thermal and zeropoint fluctuations of the atom positions, convolving with a Gaussian distribution [26,27]: where the rms x-projection of atom displacement at room temperature can be inferred directly from the dedicated experiments [28]: Finally, the potential of the entire crystal for orientation (110) is built from (75), (76) as 7 where d = 1.92Å is the interplanar distance, whereas for orientation (111), where d S = 0.78Å is the shorter interplanar distance (cf. Secs. IV B 1, IV B 2), and ⌊n/2⌋ is the (lower) integer part of n/2. One of the most commonly used parametrizations for the atomic potential screening function is that proposed long ago by Molière [29] as an analytic interpolation to the Thomas-Fermi approximation: (80) 7 Albeit at evaluation of the VR angle the integration in (9), (10) and in the subsequent equations is carried out over a single period of the crystal, in (78) it is necessary to take a sufficiently large number of planes in order to ensure the symmetry of the potential wrt any of its minima and maxima.
implies that g ′′ (r) along with the product rn e (r), where n e (r) is everywhere finite, must vanish 8 at r = 0: In contrast, Molière parametrization gives g ′′ (0) = a −2 T F α k β 2 k > 0, since all α k > 0. There is experimental evidence (see, e.g., [30]) that Molière parametrization predicts an excessive screening radius. There were proposed various heuristic recipes for reducing it [27,31], for instance [32,33], (84) The most accurate approach for calculation of atomic structure is relativistic Hartree-Fock (HF). Two commonly used parameterizations of HF calculations of atomic structure were proposed by Doyle and Turner 8 For the static continuous potential of a single plane, according to Eq. (75), that implies V ′′′ P S (0) = 0.
[34]. 9 One of them is for the electron scattering form factor with a k , b k given in Table 4 of [34], which is convenient for deriving the mean electrostatic potential: It satisfies condition (83), but unfortunately, is plagued by a more severe intrinsic flaw: g(0) = 0 instead of 1 (see Fig. 4, dotted curve), and g ′′′ (0) = 0 instead of value 4π Z n e (0) following from Poisson equation (82). This parametrization is nonetheless often used for channelingrelated problems, including VR [1,10,37].
Another Doyle-Turner parametrization is for X-ray scattering factor: where a k , b k , and c = Z − 4 k=1 a k > 0 are given in Table 3 of [34] or Table 6.1.1.4 of [38]. Constant c ∼ 1 is arguably non-physical, since it contradicts the correct asymptotic behaviour f X (s) ≃ s→∞ |n ′ e (0)| 2 5 π 3 s 4 , but parametrization (87) is intended to hold only for s ≤ 2Å −1 . In its original form, extended to arbitrary s, where erfc is the complementary error function, it was used for description of VR in [13,19]. It obeys requirement (83), but predicts g(0) = 1 − c/Z < 1 (i.e., c out of Z electrons sitting exactly at the nucleus), 10 and vanishing initial fourth derivative: g (iv) (0) = 0, instead of g (iv) (0) = 4π Z n ′ e (0) following from Poisson equation (82). A completely self-consistent parametrization of HF data, was proposed in [39]. It consists of a sum of exponentials, like (80) (permitting some of parameters a k to be negative), but is formulated for the electron density rather than the screening function. The latter derives as with and the thermal-averaged potential of a single atomic plane reads Although screening functions of HF potentials described above are markedly different at small r (see Fig. 4), the corresponding continuous planar potentials obtained after integration over longitudinal coordinates (75) and thermal averaging (76) are visually indistinguishable. 11 In Figs. 5a,b they correspond to the same solid curves. The Molière continuous potential, though, deviates from HF markedly. The reason is that Molière screening function deviates from HF at large r, whereas small-r variations are less important, being smeared out by averaging procedures, anyway. However, Molière potential with modified screening radius (84), even though fortuitously, is very close to HF, as well. In this sense, it can be regarded as competitive.

C. Comparison with experiments
The theoretical predictions for the mean VR angle computed by Eq. (12) or (13) with Molière and HF potentials, as well as NLO predictions for HF potential computed by Eqs. (28), (29), are compared with the experimental data in Figs. 6-9. Given that room-temperature planar potentials for all types of HF screening functions are very close (Fig. 5), it is natural that their predictions for the mean VR angle are indistinguishable in the figures. It can be seen from the plots that the HF potential predictions are generally the most accurate, whereas the original Molière potential, as was previously pointed out in [19], tends to give excessive predictions. But its relative deviation is twice smaller than that for the potential well depths in Fig. 5, which is natural in view of the square root dependence on the potential in Eqs. (12), (13). So, the sensitivity of χ to the interplanar potential shape can be regarded as moderate. The predictions of the Molière potential with modified screening radius (84), in fact, virtually fall on top of the HF curves, too.
The curve slopes depend on the atomic potential weakly, so, the main differences in their values stem from   Points: blue, [5]; green, [6]; red, [9]. Curves: the same as in Fig. 6.
the differences in χ 0 . The latter, according to Eqs.  Table I. This table allows one to assess subtle differences between the predictions of different potentials.
It should be cautioned that the pure continuous potential deflection theory is valid under condition where the parameter in the rhs, quantifying the onset of  (111) at R = 15 cm. Points, experimental data [10]. Curves: the same as in Fig. 6. A break in the curves around E = 18 GeV occurs because the minor potential bump rises above the major barrier due to the large tilt, making the potential double-well instead of single-well. At that, one has to switch from Eq. (12) to Eq. (13). Points: blue, measurement [7] for E = 150 GeV and R = 22.79 m; green, measurement [11] for E = 120 GeV and R = 2.71 m. Curves: the same as in Fig. 6.
incoherent multiple scattering effects, 12 may be evaluated 12 I.e., breakdown of condition where the lhs is the typical angle of atomic plane crossing by a VR particle near the radial reflection point, and the rhs is the rms incoherent multiple scattering angle acquired in the intrinsic VR region of the length ∆z defined by Eq. (7). Violation of condition (94) will lead to impossibility for the particle to pass the entire sequence of potential wells, and hence to a decrease of the reflection angle.
by the engineering formula (64) of [17]: At R −1 R −1 mult , the experimental points in Figs. 6-7 subside below the theoretical prediction. Quantitative explanation of this effect is challenging for the analytic theory.
The logarithmically modified NLO approximation developed in Sec. III C (blue dot-dashed curves in Figs. [6][7][8][9] proves to work nicely for negatively charged particles (Figs. 8-9). For positively charged particles, its application leads to somewhat worse, slightly excessive predictions (compared with the exact result) in the intermediate region. This resembles the situation with a square-well, rather than parabolic potential (see the end of Sec. IV A), although visually interplanar potentials in Fig. 5 look close to parabolic. That occurs because, as was mentioned in Sec. III C, the thermally smeared potential around the atomic planes is not described by a reverted parabola at a long extent, and at intermediate R c /R, radial reflection points belong to a non-parabolic region.
If instead one employs the NLO formula for the static potential [Eq. (25)], which displayed in Figs. 6, 7 by dotted curves, it proves to be closer to the exact result at moderate R c /R, while being somewhat worse at R c /R → 0, which is not surprising since this approximation is not truly asymptotic. But physically, at very small R c /R, as was mentioned above, the pure potential description of deflection is invalid, anyway. Therefore, in practice, the use of the static-potential NLO approximation may have certain advantages. Finally, it attracts attention that experimental data for negatively charged particles, at least for orientation (111) (see Fig. 8), overshoot the theoretical curves, not being reproduced by any of the most commonly used potentials sufficiently well (as was mentioned above, the discrepancies are primarily in χ 0 ). Furthermore, the prediction by Molière potential appears to be closer to those data than HF, although in other cases it was vice versa. In this regard, one can consider two basic possibilities.
The first one is to recall that the experimental procedure for determination of χ consists in fitting the measured angular distribution of transmitted particles by a superposition of Gaussian-related functions, each of which is associated with a certain particle fraction (channeled, dechanneled or VR). But for negatively charged VR particles, the angular distribution has a rather long 'tail' towards the crystal bending due to "orbiting" of particles at the tops of broad interplanar potential maxima [15]. Such particles can be misidentified as not belonging to the VR fraction (instead being included to the quickly dechanneled and multiply scattered fraction), and their subtraction increases the mean VR angle. 13 Secondly, it should be borne in mind that even though the noninteracting-atom approximation is rather good for third-row chemical elements (much better than for second-row elements [42,43]), covalent binding effects can nonetheless somewhat alter the interplanar potential [40][41][42], and hence have some impact on χ 0 . Recalling also (see Fig. 4) that the observed difference of the continuous planar Molière potential from planar HF potential arose due to a difference between their atomic screening functions at large r (∼ 1Å), where the applicability of the noninteracting-atom approximation generally breaks down, it is not excluded that real deformations of the valence electron density may give rise to observable differences in continuous potentials, as well.
If one considers a possibility that the good agreement for χ 0 for positively charged particles could be partly fortuitous, and looks for a modification of the interplanar potential such that the value of χ 0 is unaffected for positively charged particles, but increases for negatively charged ones, qualitatively, it may be instructive to refer to the double-parabolic potential model of Sec. IV B 3 [for simplicity, for the case of orientation (110)]. There, χ 0 for positively charged particles was expressed by an explicit function χ 0 ∝ V 0 (1 − δ) of the potential depth V 0 and shape parameter δ [see Eq. (64)]. If the latter is increased so that product V 0 (1 − δ) remains constant, for positively charged particles χ 0 will not change, whereas for negatively charged ones, it will increase by factor 1+δ [see Eq. (67)]. Physically, the increase of δ corresponds to an increase of the potential curvature midway the planes, which from the viewpoint of the Poisson equation does not contradict to formation of covalent bonds in that region. But to achieve quantitative agreement with the data, one needs sizeable δ ∼ 0.1, whereas the fitted value 13 It may be mentioned yet that experiment [10] was performed with a rather thin crystal (L = 60 µm), whereas herein we presume validity of the thick-crystal limit. Nonetheless, condition (8) for it is fulfilled. Besides that, boundary effects are not expected to be strictly positive, and thus cannot explain the stable excess of the data over the theory predictions.
(66) is significantly lower. Proper investigation of this issue is beyond the scope of the present paper.

VII. SUMMARY AND OUTLOOK
The developed approach for evaluation of the mean volume reflection angle solves the problem of summation over the crystal periods exactly. Apart from the apparent convenience for making numerical predictions, it also permits a number of qualitative inferences. In particular, it reveals a duality between VR and channeled particles, in the sense that only under-barrier regions of the effective potential contribute to the mean VR angle. The limiting (at 4R c /R ≪ 1) value of the latter proves to be related with a net critical channeling angle for a straight crystal (its local generalization averaged over the interplanar interval).
It is also noteworthy that in practice, for description of E-and R-dependencies of χ, it may be sufficient to use NLO expansion in the small ratio R c /R (Sec. III). That demands the knowledge of only a few empirical constants, calculable numerically based on a realistic interplanar potential for a straight crystal. NLO formulas for orientation (110) and for parabolic interplanar potential, were formerly implemented in the Monte-Carlo code FLUKA [46]; they can be improved based on the present development.
Cases of medium R/4R c and R → R c were also studied. At R 4R c , some guidance can be provided by parabolic potential models of Sec. IV, offering simple closed-form expressions for χ. The most accurate predictions, though, are obtained with HF potentials.
The comparison between the theory and results of experiment [10] on VR of negatively charged particles and crystal orientation (111) displays yet some discrepancy (Sec. VI C), the interpretation of which requires further investigation. More data on VR of negatively charged particles, for all crystal orientations are desirable.
It is worth recapitulating that in the present paper, VR was everywhere treated as a purely elastic deflection process. In reality, it is accompanied by weakly inelastic and quasi-elastic scattering processes, which need to be understood, as well. In particular, they give rise to a non-negligible volume capture probability, due to which the VR efficiency even at high energy can amount to less than 90%.
Another kind of inelastic processes is emission of electromagnetic radiation, which can be intense for ultrarelativistic electrons and positrons passing through crystals [8,11,25]. For radiation at VR, the deviation of the spectrum from the kinematically formed coherent bremsstrahlung in a bent crystal (generated away from the intrinsic VR region) is concentrated in the low-ω region [44]. The infrared limit ω → 0 itself is rather trivial, corresponding to factorization of the radiation and scattering probabilities [44]. Less trivial is the slope of the radiation spectrum at the origin. For the latter quan-tity, there was derived a generic expression depending on the detail of the particle trajectory inside the target [45]. To describe that behavior, however, one needs to go beyond the present theory, which deals only with the final deflection angle.
Appendix A: Derivation of formula (10) In this appendix we derive formula (10) central for the present paper. If one is interested in the mean VR angle alone, rather than in the detail of the angular distribution, it can be evaluated by interchanging the order of integrations in the double integral arising when Eq. (5) is inserted into Eq. (4). The dependence of the integrand on E ⊥ is very simple and should thus be manageable exactly.
To interchange the order of integrations, one must first take into account that the lower limit r min (E ⊥ ) of r integrations is E ⊥ dependent. Since r min (E ⊥ ) by definition is a maximal r such that V eff (r) = E ⊥ , condition r > r min (E ⊥ ) equivalently expresses as an r-dependent lower bound for E ⊥ : The absolute minimum of r is subsequently inferred by noting that the dependence (A1) is monotonic, whence the minimal r corresponds to the maximal E ⊥ , i.e., to const + ∆E ⊥ . Without the loss of generality, the origin of coordinate r can be chosen so that r min (const + ∆E ⊥ ) = 0. Therewith, The corresponding integration domain in the {r, E ⊥ } plane is shown by a horizontal band in Fig. 10.
We are now in a position to make use of the periodicity of V (r). The integrand of the double integral in (A2) is periodic both in r and E ⊥ at constant E ⊥ + E R r. If we split the integral with respect to r into a sum of integrals over full periods: shifts in E ⊥ , whereby we obtain a sum of integrals over E ⊥ : One observes that they seamlessly combine into a single integral with respect to E ⊥ , with the upper limit ∆E ⊥ r 0 /d = Er 0 /R: Graphically this transformation is illustrated in Fig. 10. The E ⊥ -integration then performs exactly in a trivial manner, and the limit r 0 → ∞ is evaluated to give Eq. (10).
Appendix C: Derivation of formula (64) Here we derive an approximation for the mean VR angle in a single-well potential described by two joint parabolas, Eq. (55), suitable under conditions (62) and d 1 ≪ d 2 . Our starting point is Eq. (12).
When the centrifugal potential − Er R is added to the crystal potential (58), under condition R ≫ R c , the location of the maximum of the tilted effective potential shifts from d only slightly: r m = d − The effective potential value in this maximum equals V eff (r m ) = V 3 (r m ) − Erm R . Exact evaluation of the integral (12) for potential (55) then gives or, evaluating the latter simple integrals, Invoking condition (59), wherewith Ed1 4V1 = Ed2 4V2 = R c , expression (C2) simplifies to Notably, all the E/R dependence here is given by the same factor 1 − Rc R 2 as for the purely parabolic well, cf. Eq. (33).
In the limit R/R c → ∞, however, expression (C3) does not exactly tend to χ 0 , because it was derived presuming condition (62), where the right-hand side is large but finite. Its precise limiting value can be obtained by returning to formula (14): or, with the use of condition (59), But as long as V 1 /V 2 0.1, the difference That approximation can already be used for arbitrarily large R.