1 Introduction

There exists since long considerable interest in the channeling process of ultra-relativistic electrons and positrons at planes of a single crystal. Of particular interest is the emission of undulator-like radiation in periodically bent crystals aiming in the construction of compact radiation sources in the MeV range and beyond, see, e.g. Korol et al. [1]. However, collisions of the leptons with atoms and their electrons comprising the crystallographic planes lead to a de-channeling process being in conflict with the envisaged application. It is, therefore, of utmost importance to understand the channeling as well as the de- and also re-channeling processes in detail, in particular also for electrons.

While “de-channeling” is a well-defined technical term, the de-channeling length is not. Indeed, a particle entering a straight channel needs some time, or a certain length interval, to equilibrate. Therefore, a constant de-channeling rate, leading to an exponential decrease of the channel occupation does generally not exist. Reliable experimental information on the rate distribution as function of the penetration depth of the particle in the crystal is rare, or does not exist at all. However, such information is required to determine the channeling efficiency, i.e., that fraction of particles which follow the exponential decay law. This paper intends to contribute to this task via Monte Carlo simulations. This way, also the re-channeling process has been studied, including the time when re-channeling happens, the number of planes the electron passes until it finally re-channels, the probability distribution, and the efficiency of the re-channeling process.

For many years channeling of ultra-relativistic particles in single crystals was described in the continuum potential picture introduced by Lindhard [2]. In particular, the Fokker–Planck equation with which the de-channeling process can be calculated is based on it, see, e.g. the papers of Backe et al. [3, 4] and references cited therein. Objections against the reliability of the solutions of the Fokker–Planck equation have been formulated by Tikhomirov [5, ch. 3.4] but it has not been investigated whether its predictive power can be improved, or not. This paper will contribute also to this subject, in particular it focuses on a better approximation of the Kitagawa–Ohtsuki drift and diffusion coefficients [6] still in use.

If experimental results do not exist, insight into the addressed issues can be obtained by Monte Carlo simulation calculations. There exist detailed considerations with sophisticated codes, see, e.g. Pavlov et al. [7] and references cited therein, and the continuum potential picture seems to be outdated. However, a detailed comparison between the former and the latter with a discussion of advantages and disadvantages does not exist. This issue will be addressed in passing.

The continuum potential picture has the advantage that channeling can be classically understood for ultra-relativistic particles in a rather intuitive manner. In particular, for ultra-relativistic leptons the dynamics can be treated classically. Deficiencies which were already pointed out by Lindhard, namely that the potential has fluctuations in the Angstrom scale, should have only little impact. A lepton moving in forward direction at angles in the order of sub milli-radians feels during its passage the potential of thousands of single atoms. It is therefore expected that the potential can be described rather precisely in such a collective approach. Moreover, it can be considered as decoupled from hard scattering events at atoms and bound electrons since those interactions are rare.

Channeling in the continuum potential picture was treated by a number of authors, for an overview see, e.g., Korol et al. [8, Chapter 2]. In particular, the DYNECHARM++ code of Bagli and Guidi [9], and CRYSTALRAD of Sytov et al. [10] should be mentioned. In this work yet another code was developed which is restricted to the above-mentioned problems, however, with the possibility to investigate the sensitivity on special interactions like low energy plasmon excitations.

The paper is organized as follows. Sections 25 summarize the theoretical background. In Sect. 2 basic definitions and the particle dynamics are described. In Sect. 3 the Doyle–Turner representation of the electron scattering factors at screened atoms is introduced with a modification into the Molière approach. The latter has the advantage that it avoids a cutoff at large momentum transfers. On the basis of this data set the electric potential, the electron density, and collision cross-sections are calculated. The latter finally leads to the scattering distribution needed for the Monte Carlo simulations. In Sect. 4 and in the appendix A the electron-electron interaction is treated. At energy losses larger than about 30 keV the Møller cross section has been applied. For less than 30 keV models of Ashley [11] and Fernández-Varea [12] were utilized leading to the double differential cross-sections as function of both, the energy and momentum transfer to bound electrons from which the electron-electron scattering distribution was derived. In Sect. 6 preparatory calculations are described which include the test of the model using the example of an amorphous carbon foil, and the calculation of the initial probability distribution. In Sect. 7 results of the Monte Carlo simulations for channeling and re-channeling of 855 MeV electrons in straight and bent (110) planes of diamond are presented, and for the latter compared with published results obtained with the MBN explorer package. In Sect. 8 drift and diffusion coefficients, which enter the Fokker–Planck equation, were calculated for both, scattering at atoms and separately at electrons, and the impact on its predicting power is discussed. The paper closes with a conclusion in Sect. 9.

2 Particle dynamics

For a plane crystal a laboratory coordinate system (xyz) has been chosen with the z-axis located in the crystallographic plane, the horizontal x-axis perpendicular to it, and the vertical y-axis such that a right-handed coordinate system results. The planes of a bent crystal are assumed to be circularly shaped with a constant bending radius R, and a concave shape with respect to the z-axis. The entrance angle into the crystal \(\psi _0\) is the projected angle onto the (x,z) plane measured with respect to the z axis. The rotation of the crystal around the y-axis has a negative sign for clockwise, and positive one for counter clockwise rotation. Analytically the following expressions were used which were derived from the equation of a circle with its centre at \(x_0 =- R \cos \psi _0 \) and \(z_0 = R \sin \psi _0 \):

$$\begin{aligned}{} & {} x = R \cdot \left( \sqrt{1-(\sin \psi _0-z/R)^2}-\cos \psi _0~\right) \end{aligned}$$
(1)
$$\begin{aligned}{} & {} \frac{\textrm{d}x}{dz} = \frac{\sin \psi _0-z/R}{\sqrt{(1-(\sin \psi _0-z/R)^2}} \end{aligned}$$
(2)
$$\begin{aligned}{} & {} \frac{d^2x}{dz^2} \cong -\frac{1}{R}{~~~\text{ for }~~\psi _0\ll 1}. \end{aligned}$$
(3)

In the simulations it has been assumed that the particle moves freely a certain distance \(\varDelta z_k\), for definition see Eq. (4) below, which maybe randomly interrupted by hard collisions with lattice atoms and/or electrons. In the free movement the total energy is conserved while in a collision a sudden transverse energy transfer happens without a change of the x coordinate. The calculations were performed in the \((x, E_\bot )\) phase space in a grid with step sizes \(\varDelta x \simeq \) 0.01 Å  and exact transverse energy \(E_\bot \) in the potential well. For the movement of a particle at position \(x_k\) with transverse energy \(E_{\bot ,k}\) to \(x_{k+1}~=~x_k + \varDelta x_k\) a time \(\varDelta t_k\) is required. In the following the time is always multiplied by the particles velocity \(v\simeq c\), and one obtains for the step size in z direction

$$\begin{aligned} \varDelta z_k=\varDelta t_k~v=\int \limits _{x_k}^{x_k+\varDelta x_k}\frac{\sqrt{p v/2}}{\sqrt{E_{\bot ,k}-U(x)}}\textrm{d}x. \end{aligned}$$
(4)

Here \(p v = \beta \gamma m_e c^2\) is the momentum of the particle, with \(\gamma =1/\sqrt{1-\beta ^2}\) the Lorentz factor, \(\beta =v/c\) and \(m_e\) the electron rest mass. Despite the rather small integration interval \(\varDelta x_k\), the integral representation must be chosen to obtain reliable \(\varDelta z_k\), in particular at the turning points in the potential wall U(x) since the denominator tends there to zero. In that case \(\varDelta x_k\) reduces to \(\varDelta x_k = \mid U^{-1}(E_{\bot ,k})-x_k\mid \) and \(\varDelta z_k\) doubles. Eq. (4) is written for a bent crystal for which

$$\begin{aligned} U(x)=u(x)-(pv/R)x, \end{aligned}$$
(5)

with u(x) the potential of the plane crystal. However, this approach only holds for a constant bending radius R, i.e., not if the bending radius is a function of z.

At the distance \(\varDelta z_k\) a mean number of collisions happen, given by the total cross-section, from which according to the Poisson-statistics integer numbers are randomly generated. Thereafter, the individual scattering angle is randomly simulated for each scattering event with a procedure described below from which the transverse energy change is calculated.

Fig. 1
figure 1

Electron scattering factors for atomic carbon. The black full curve \(f_e^{M}\) represents the Molière approximation [13] according to Eq. (8), the red one \(f_e^{DT}\) that one of Doyle and Turner [14] in the 6-parameter representation of Eq. (9) with values taken from Ref. [15]

3 Interactions with atoms

For the calculation of the channeling potential, as well as for the also needed angular distribution of the scattered electrons, the electron scattering factors at screened atoms are required. There exist essentially two approaches to access them, that of Molière [13] and that of Doyle and Turner [14]. Using the parameters

$$\begin{aligned} \alpha ^M= & {} \{0.1,~~0.55,~~0.35\}\end{aligned}$$
(6)
$$\begin{aligned} \beta ^M= & {} \{6.0,~~~1.2,~~~~~ 0.3\}/a_{TF} \end{aligned}$$
(7)

with which Molière adapted the Thomas-Fermi function, a scattering factor

$$\begin{aligned} f_e^{M} (s)=\frac{4\pi Z \alpha \hbar c}{a^3} \sum _{i=1}^{3}\frac{\alpha _i^M}{(4\pi s)^2+(\beta _i^M)^2} \end{aligned}$$
(8)

results. It is compared in Fig. 1 with the Doyle–Turner representation

$$\begin{aligned} f_e^{DT} (s)=\frac{2\pi a_0 \alpha \hbar c}{a^3} \sum _{i=1}^{6}a_i \exp (-b_i s^2) \end{aligned}$$
(9)

using six pairs of parameters (\(a_i,b_i\)) which were quoted by Chouffani and Überall [15]. Here are \(a_{TF}=0.8853 ~a_0 ~Z^{-1/3}\) the Thomas-Fermi screening factor, \(a_0\) the Bohr radius, Z=6 for carbon, and a = 3.5668 Å the lattice constant of diamond [16]. The quantity s is related to the momentum transfer by \(q=2p v/(\hbar c) \sin (\vartheta /2) = 4\pi s\), with \(\vartheta \) the scattering angle. Significant deviations up to a factor of 2.6 can be observed for \(s<0.2 \)Å\(^{-1}\) while for larger values they are less than 8%. Nevertheless, all calculations in this paper have been performed in the Molière approximation, however, with the modified parameter set

$$\begin{aligned} \alpha= & {} \{0.09567,~2.9992\cdot 10^{-5},~~0.90430\}\nonumber \\ \beta= & {} \{9.61747, ~1.2,\qquad \qquad \quad 0.73393\}/a_{TF} \end{aligned}$$
(10)

which approximates the Doyle–Turner representation to better than 8% in the full range of \(0 \le s/\)Å \(\le 6\). The reason is that \(f_e^M(s)\) asymptotes to a Lorentzian rather than to a Gaussian, avoiding the unrealistic cut-off for large momentum transfers.

3.1 Calculation of the continuum potential

The potential u(x) has been calculated according to chapter 9.1 of the textbook of Baier et al. [17] with the Fourier-expansion coefficients

$$\begin{aligned}{} & {} G(q)=\frac{4\pi Z \alpha \hbar c}{a^3} \exp \left( -\frac{u_1^2 q^2}{2}\right) S(q) \sum _{i=1}^{3}\frac{\alpha _i}{q^2+\beta _i^2} \end{aligned}$$
(11)

and S(q) the structure factor. The one-dimensional thermal vibration amplitude is \(u_1\) = 0.04226 Å [18]. The potential u(x) is shown in Fig. 2a, and in panel (b) the electron density \(n_{el}(x)\) across the channel. The latter has been derived from the potential u(x) with the aid of the one-dimensional Poisson equation.

Fig. 2
figure 2

a Potential of the (110) plane for a plane diamond crystal as function of the inter-planar distance coordinate x. The potential minimum has been normalized to zero, i.e. the continuum border is located at \(E_{\bot ,0}\) = \(u_0\) = 23.957 eV. b Electron density \(n_e(x)\). The interplanar distance is \(d_p~=~1.261~\) Å  and the mean electron density for diamond \((1/d_p)\int _{-d_p/2}^{d_p/2} n_e(x) \textrm{d}x = Z\cdot (8/a^3)\) = 1.0578/Å\(^3\)

3.2 Scattering distribution for atoms

Collision with atoms will be described in the screened potential with the modified Molière parameters of Eq. (10). The differential cross-section reads

$$\begin{aligned}{} & {} \frac{\textrm{d}\sigma ^{(at)}}{\textrm{d}\varOmega }(\vartheta )\nonumber \\ {}{} & {} \quad =4(Z\alpha )^2 \left( \frac{\hbar c}{pc}\right) ^2 \left( \sum _{i=1}^{3}\frac{\alpha _i}{4 \sin ^2(\vartheta /2)+(\beta _i\cdot \hbar c/pc)^2}\right) ^2. \nonumber \\ \end{aligned}$$
(12)

Eq. (12) can be derived from the Mott cross-section by replacing \(1/q^2\) by

$$\begin{aligned} \frac{1-F(q)}{q^2}=\sum _{i=1}^{3}\frac{\alpha _i}{q^2+\beta _i^2}, \end{aligned}$$
(13)

with F(q) the atomic form-factor. At planar channeling we are interested in the differential scattering cross-section with respect to the x coordinate. It is obtained from Eq. (12) for \(\vartheta \ll 1\), and \(\vartheta ^2=\vartheta _x^2+\vartheta _y^2\) by integration over the \(\vartheta _y\) coordinate and reads approximately for \(\vartheta _x\ll 1\)

$$\begin{aligned}{} & {} \frac{\textrm{d}\sigma ^{(at)}}{\textrm{d}\vartheta _x}(\vartheta _x)\cong 4(Z\alpha )^2 \left( \frac{\hbar c}{pc}\right) ^2\times \nonumber \\{} & {} \qquad \quad \int _{-\infty }^\infty \left( \sum _{i=1}^{3}\frac{\alpha _i}{\vartheta _x^2+{\vartheta _y'}^2+ (\beta _i\cdot \hbar c/pc)^2}\right) ^2 d{\vartheta _y}'.\nonumber \\ \end{aligned}$$
(14)

From the total cross-section

$$\begin{aligned} \sigma _{tot}^{(at)}=\int _{0}^\pi \frac{\textrm{d}\sigma ^{(at)}}{\textrm{d}\varOmega }(\vartheta )2\pi \sin \vartheta d \vartheta \cong \int _{-\infty }^\infty \frac{\textrm{d}\sigma ^{(at)}}{\textrm{d}\vartheta _x}(\vartheta _x) \textrm{d}\vartheta _x \end{aligned}$$
(15)

one obtains the mean number of collisions in an interval \(\varDelta z_k=\varDelta t_k~v\) as

$$\begin{aligned} \varDelta m_k=\frac{8}{a^3}\varDelta z_k~\cdot \sigma _{tot}^{(at)} \frac{d_p}{\sqrt{2\pi }u_1} \exp (-x_k^2/2u_1^2). \end{aligned}$$
(16)

Here a = 3.5668 Å  is the lattice constant of diamond [16], \(8/a^3\) the atomic number density, and \(d_p=a/2\sqrt{2}\) = 1.261 Å  the inter-planar distance. Due to the distribution of the atoms comprising the planes by thermal vibrations, the mean number of collisions is also a function of the inter-planar coordinate \(x_k\). The probability distribution for m collisions at a mean \(\varDelta m_k\) is given by the Poisson-distribution function

$$\begin{aligned} P_m(k)=\frac{(\varDelta m_k)^m}{m!}\exp (-\varDelta m_k). \end{aligned}$$
(17)

The underlying assumption for applying this Eq. (17) is statistical independence of the collisions which should be fulfilled if the particle direction does not coincide with a crystal axis. In the simulation \(M_k = \text{ Random}_{m}[P_{m}(k)]\) is randomly generated. For \(M_k\) = 0 the change of the scattering angle is \(\varDelta \vartheta _{k}^{(at)}\) = 0, otherwise

$$\begin{aligned} \varDelta \vartheta _{k}^{(at)} = \sum _{i=1}^{M_k} \text{ Random}_{i,\theta _{x}}[P^{(at)}(\theta _{x})], \end{aligned}$$
(18)

is obtained by a repeated generation of random numbers from the normalized scattering distribution function \(P^{(at)}(\theta _{x})\). The latter has been derived from the angular part of Eq. (14) after normalization with the total cross-section Eq. (15). It is a rather involved function of the parameters \(\alpha _i\) and \(\beta _i\). Numerically one obtains

$$\begin{aligned} P^{(at)}(\theta _x)= & {} \left( \frac{1.68969 \cdot 10^{-8}+\theta _x^2}{(4.3114\cdot 10^{-11} + \theta _x^2)^{3/2}} \right. \nonumber \\{} & {} -\frac{3.56831 \cdot 10^{-12} +0.03096\cdot \theta _x^2}{(1.1526\cdot 10^{-10} + \theta _x^2)^{3/2}} \nonumber \\{} & {} \left. -\frac{6.98552 \cdot 10^{-9} + 0.96904\cdot \theta _x^2}{(7.4034\cdot 10^{-9} + \theta _x^2)^{3/2}}\right) \Big /786.89.\nonumber \\ \end{aligned}$$
(19)

The distribution function is shown in Fig. 3, full curve. The total scattering cross-section according to Eq. (15) is \(\sigma _{tot}^{(at)}\) = 24.45\(\cdot 10^{-4}~\)Å\(^2\). A mean transverse energy gain \(<\varDelta E_\bot /\varDelta z>_{at}\) = \(8/a^3 \sigma _{tot}^{(at)} pv/2 <\theta _x^2>_{at}\) = 0.709 eV/\(\mu \)m is obtained for a somehow arbitrarily chosen cut-off angle of 0.02 rad. The mean number of collisions is \(8/a^3 \sigma _{tot}^{(at)}\) = 4.31/\(\mu \)m.

Fig. 3
figure 3

Normalized atomic (full) and electronic (dot-dashed) scattering distributions \(P^{(at)}(\theta _x)\) and \(P^{(el)}(\theta _x)\), respectively, as function of the scattering angle \(\theta _x\) for 855 MeV electrons on carbon atoms. The electronic distribution is very narrow extending vertically up to 3.73\(\cdot 10^6\). For a logarithmic representation of the electronic distribution see Fig. 16. The FWHM amount to 10.08 \(\mu \)rad (atomic), and 0.0752 \(\mu \)rad (electronic). Both distributions have long tails taken into account in the numerical simulation up to ± 0.02 rad (atomic) and ± 0.0345 rad (electronic). The root mean squared scattering angles are \(<\theta _x^2>^{1/2}\) = 19.6 \(\mu \)rad (atomic), and 7.33 \(\mu \)rad (electronic)

4 Electron-electron interactions

4.1 General remarks

The mean transverse energy gain \(\overline{\varDelta E_\bot ^{el}}\) due to the scattering of the beam particle at atomic electrons in an interval \(\varDelta z=\varDelta t~v\) can approximately be related to the mean electronic energy loss by the relation [19]

$$\begin{aligned} \frac{\overline{\varDelta E_\bot ^{(el)}}}{\varDelta z}=\frac{\alpha _P}{2\gamma }\frac{\overline{\varDelta E^{(el)}}}{\varDelta z} \end{aligned}$$
(20)

with \(\alpha _P~\simeq ~0.5\). The mean energy loss \(\overline{\varDelta E^{(el)}/\varDelta z}\) can either be taken from tables, e.g. [20], or calculated for higher energies than 1 GeV with the equation [19]

$$\begin{aligned} \frac{\overline{\varDelta E^{(el)}}}{\varDelta z} = \frac{2\pi (\hbar c \alpha )^2}{m_e c^2 \beta ^2}\frac{8 Z}{a^3}L_e, \end{aligned}$$
(21)

with the Sternheimer correction factor [21, Eq. (52), Eq. (2), Table 2]

$$\begin{aligned} L_e=\ln \left( \frac{\gamma ^2 2 m_e c^2 pv/2}{I^2}\right) -2\beta ^2-4.606\cdot \lg (\gamma )+4.636.\nonumber \\ \end{aligned}$$
(22)

The quantity I is the ionization potential for which I = 89.4 eV has been chosen [22], \(\alpha \) the fine-structure constant, and \(\beta ~=~v/c\). With \(\overline{\varDelta E^{(el)}/\varDelta z}\) = 737.82 eV/\(\mu \)m, from Eq. (21), one obtains from Eq. (20) \(\overline{\varDelta E_\bot ^{(el)}/\varDelta z}\) = 0.110 eV/\(\mu \)m.

Unfortunately, this approach does not lead to a viable simulation scheme. First of all, Eq. (20) describes a mean value, however, required is the underlying scattering distribution as function of the scattering angle. In addition, \(\alpha _P~\simeq ~0.5\) is a rather crude estimate.

There are different possibilities to proceed. The simplest one is to neglect scattering on electrons, leading to an underestimate of the transverse energy transfer. Another at the first glance also simple one would be to treat the electrons as a free gas and apply Møller scattering. However, a proper low energy cutoff is required since the Møller cross-section diverges at zero energy transfer. The third possibility applied here requires model assumptions on the generalized oscillator strength of the atom, as described below and in more detail in appendix A. In the latter, also the implications of the free electron gas approach are discussed in more detail.

4.2 Scattering distribution for electrons

From Eq. (51) of appendix A one obtains for the angular differential cross-section

$$\begin{aligned} \frac{\textrm{d}\sigma ^{(el)}}{\textrm{d}\varOmega }(\theta )= & {} \int \limits _{W_{min}(\theta )}^{W_{m}}dW\frac{d^2 \sigma ^{(el)}}{\textrm{d}W\textrm{d}\varOmega }(\theta ,W)~~\varTheta (\theta _{m}-\theta )+ \nonumber \\{} & {} \frac{\textrm{d}\sigma _{Mo}}{\textrm{d}\varOmega }(\theta )~~\varTheta (\theta -\theta _{m}). \end{aligned}$$
(23)

with the matching energy \(W_{m}\) = 27 keV, and the matching angle \(\theta _{m} = \theta _{max}\)(\(W_m\)) = 194.2 \(\mu \)rad from righthand side of Eq. (53). For \(W > W_{m}\) all electrons of the carbon atom are regarded as quasi-free, and the Møller cross-section in the laboratory system, Eqs. (1) and (2), of Ref. [23] can be applied (read in Eq. (1) for the prefactor of the last term \((\gamma -1)^2/\gamma ^2\) instead of \((\gamma -1)^2/\gamma \)). It should be mentioned, that at \(W_{m}\) the cross-sections, first and second term of Eq. (23), match with an accuracy of 1.9 %. This is a remarkable result in view of the fact that both cross-sections were calculated by complete different approaches.

The total cross-section is the integral

$$\begin{aligned} \sigma _{tot}^{(el)} = \int \limits _{\vartheta =0}^{\varTheta _{max}} \textrm{d}\vartheta ~2\pi \sin \vartheta \frac{\textrm{d}\sigma ^{(el)}}{\textrm{d}\varOmega }(\vartheta ). \end{aligned}$$
(24)

For \(\varTheta _{max} = \arcsin \sqrt{2/(\gamma +3)}\)= 0.03454 rad has been chosen, corresponding to the maximum energy loss of \((\gamma -1)m_e c^2/2\) = 427.5 MeV. The total cross-section for scattering at a single electron is \(\sigma _{tot}^{(el)}\) = 5.576 \(\cdot 10^{-4}~\)Å\(^2\). This result agrees rather precisely with the integral \(\int (\textrm{d}\sigma ^{(el)}/dW)~dW = 5.577 \cdot 10^{-4}~\)Å\(^2\) of Eq. (49) in appendix A.

Again, at planar channeling we are interested in the differential scattering cross-section with respect to the x coordinate. It is obtained by integration of Eq. (23) over the \(\vartheta _y\) coordinate and reads

$$\begin{aligned} \frac{\textrm{d}\sigma ^{(el)}}{\textrm{d}\vartheta _x}(\vartheta _x) = \int \limits _{\vartheta _y=-\varTheta _{y}}^{\varTheta _{y}} \textrm{d}\vartheta _y \frac{\textrm{d}\sigma ^{(el)}}{\textrm{d}\varOmega }\left( \sqrt{\vartheta _x^2+\vartheta _y^2}\right) \end{aligned}$$
(25)

with \(\varTheta _{y}=\sqrt{\varTheta _{max}^2-\vartheta _x^2}\). The scattering distribution is the normalized angular distribution of Eq. (25)

$$\begin{aligned} P^{(el)}(\theta _x)=\frac{\textrm{d}\sigma ^{(el)}}{\textrm{d}\theta _x}(\theta _x)/\sigma _{tot}^{(el)} \end{aligned}$$
(26)

which is depicted also in Fig. 3, dot-dashed line. It can be approximated by the heuristic functionFootnote 1

$$\begin{aligned} P^{(el)}(\theta _x)=\frac{A}{|\theta _x^3 |+ c \cdot \theta _x ^2 + b\cdot |\theta _x|+ a}. \end{aligned}$$
(27)

The mean squared scattering angle

$$\begin{aligned} <\theta _x^2>_{el}~ = \int \limits _{-\varTheta _x}^{\varTheta _x} \theta _x^2 \cdot P^{(el)}(\theta _x)\textrm{d}\theta _x \end{aligned}$$
(28)

is a function of the integration limit \(\varTheta _x\). For \(\varTheta _{max}\) = 0.03454 rad \(<\theta _x^2>_{el} = 5.370 \cdot 10^{-11}\) results. For the calculation of the mean transverse energy gain, which is of relevance for channeling, one obtains \(<\varDelta E_\bot /\varDelta z>_{el}\) = \((8~Z/a^3)\) \(\sigma _{tot}^{(el)}pv/2 <\theta _x^2>_{el}\) = 0.136 eV/\(\mu \)m which would imply in Eq. (20) an improved \(\alpha _P = 0.618\). This mean transverse energy gain for electron-electron collisions is a factor of 5.2 less than the corresponding number for electron-atom collisions. The mean number of collisions per atom and unit longitudinal length interval is \((8Z/a^3)\cdot \sigma _{tot}^{(el)}\) = 5.899/\(\mu \)m,

The electron density \(n_{el}(x)\) across the channel is a function of the coordinate x. It has been derived from the potential u(x) with the aid of the one-dimensional Poisson equation. The result is shown in Fig. 2b. The mean number of collisions at the lateral position \(x_k\) in an interval \(\varDelta z_k=\varDelta t_k~c\) is

$$\begin{aligned} \varDelta n_k=n_{el}(x_k)~\sigma _{tot}^{(el)}~\varDelta z_k. \end{aligned}$$
(29)

For this equation it has been assumed that the number of collisions depends only on the electron density, i.e., a possible dependence on the transverse coordinate x of the cross-section has been neglected by the replacement with its mean \(\sigma _{tot}^{(el)}\). However, considering a certain fraction of the electrons as quasi-free “sea” electrons, say 20 %, for which the cross section can be described by Møller scattering with a low energy cutoff of 6 eV, does not change the final simulation results significantly. The probability distribution for n collisions at a mean \(\varDelta n_k\) is given by the Poisson-distribution function

$$\begin{aligned} P_n(k)=\frac{(\varDelta n_k)^n}{n!}\exp (-\varDelta n_k). \end{aligned}$$
(30)

In the simulation \(N_k = \text{ Random}_{n}[P_{n}(k)]\) is randomly generated. For \(N\ge 1\) the change of the scattering angle is

$$\begin{aligned} \varDelta \vartheta _{k}^{(el)}= \sum _{i=1}^{N_k}\text{ Random}_{i,\theta _x}[P^{(el)}(\theta _x)]. \end{aligned}$$
(31)

5 Transverse energy transfer at collisions

The Monte Carlo simulated change of the scattering angle is with Eqs. (18) and (31) \(\varDelta \vartheta _k=\varDelta \vartheta _{k}^{(at)} +\varDelta \vartheta _{k}^{(el)}\), and the new scattering angle at the step from \(x_k\rightarrow x_{k+1}=x_k+\varDelta x_k\) is given by

$$\begin{aligned} \vartheta _{k+1} =\vartheta _k+\varDelta \vartheta _{k}= \pm \sqrt{2\left( E_{\bot ,k}-U(x_{k})\right) /pv}+\varDelta \vartheta _{k} \end{aligned}$$
(32)

with \(U(x_{k})\) the potential energy at position \(x_k\). A scattering by the angle \(\varDelta \vartheta _{k}\) results in the new transverse energy

$$\begin{aligned} E_{\bot ,k+1}= & {} E_{\bot ,k}+pv/2\cdot (\varDelta \vartheta _k)^2 \nonumber \\{} & {} \pm \sqrt{2pv\cdot \left( E_{\bot ,k}-U(x_{k})\right) }~~\varDelta \vartheta _k. \end{aligned}$$
(33)

The plus sign in front of the square root holds for \(\varDelta x_k>0\), or if \(E_{\bot ,k+1}-U(x_{k+1})\) is negative, the minus sign for \(\varDelta x_k<0\). How these equations come about is explained in Appendix B. Since \(\varDelta \vartheta _k\) has both signs which are distributed with equal probability, the first term in the second row of Eq. (33) \(\varDelta E_{\perp ,k}^{(\text {drift})}\) results in a drift while the second one \(\varDelta E_{\perp ,k}^{(\text {diff})}\) in fluctuations of \(E_{\bot }\). Therefore, both terms contribute to de-channeling while re-channeling is effected only by the second term.

With these definitions instantaneous drift and diffusion coefficients may be defined as

$$\begin{aligned} \tilde{d}_1(E_{\perp ,k},~x_k)=\varDelta E_{\perp ,k}^{(\text {drift})}/\varDelta z_k=p v/2 \cdot (\varDelta \vartheta _k)^2/\varDelta z_k, \end{aligned}$$
(34)

and

$$\begin{aligned} \tilde{d}_2(E_{\perp ,k},~x_k)= & {} 1/2 \cdot (\varDelta E_{\perp ,k}^{(\text {diff})})^2/\varDelta z_k = \nonumber \\= & {} 1/2\cdot 2 p v \left( E_{\bot ,k}-U(x_{k})\right) \cdot (\varDelta \vartheta _k)^2/\varDelta z_k.\qquad \end{aligned}$$
(35)

In contrast to the coefficients \(D_1(E_\perp )\) and \(D_2(E_\perp )\) which enter the Fokker–Planck equation as mean values with respect of one oscillation period for a given \(E_\perp \), these instantaneous quantities are functions of both, the variables \(E_\perp \) and x. In addition, they have a functional dependence on the penetration depth z. We come back to this fact in chapter 8.

6 Preparatory calculations

6.1 Test of the simulation model for “amorphous diamond”

The formalism described above has been checked by a simulation of scattering of 855 MeV electrons passing an amorphous carbon foil with the density of diamond and a thickness of 120 \(\mu \)m. The mean number of collisions amounted to 517.6 and 658.3 for collisions with the atoms and electrons, respectively. The quantities \(\varDelta \vartheta _{k}^{(at)}\) and \(\varDelta \vartheta _{k_{}}^{(el)}\) were calculated with Eqs. (18) and (31), respectively. The scattering distribution as function of \(\theta _x=\varDelta \vartheta _{k}^{(at)}+\varDelta \vartheta _{k}^{(el)}\) is shown in Fig. 4. It should be mentioned that the width of the electronic distribution amounts to about 35 % of the width of the atomic distribution which cannot be neglected. The total distribution exhibits a slight asymmetry, clearly recognizable when compared with the also shown Gaussian angular distribution, which is typical for a Poisson distribution. In conclusion, the overall agreement is quite good. This fact provides some confidence in trusting the simulation calculations for the de-channeling length to be described in Sect. 7.

Fig. 4
figure 4

Simulation of the scattering distribution at the passage of 855 MeV electrons through an amorphous carbon foil with the density of diamond and a thickness of 120 \(\mu \)m. A number of 5000 tracks were simulated, error bars are of \(\sqrt{N}\) type, with N the counts in a channel. The full curve represents the projected Gaussian angular distribution according to the Particle Data Group with \(\theta _{plane}^{rms}=0.368\) mrad, which is quite similar to that one published by G.R. Lynch and O.I. Dahl [24] with 98% in the sample and \(\theta _{plane}^{rms}=0.361\) mrad

6.2 Initial occupation probability

For the initial occupation probability in the potential pocket a uniform distribution of the electron density across the transverse x coordinate, and a Gaussian scattering distribution with standard deviation \(\sigma '_{x}\) = 12.7 \(\mu \)rad for the angular divergence were assumed. The initial occupation probability can be analytically calculated according to Backe et al. [3, Eq. (15)]. The initial distribution is shown in Fig. 5 for a straight diamond crystal. Part of the distribution is located in the continuum, i.e., about 12 % occupy states with \(u>u_0\) = 23.96 eV.

Fig. 5
figure 5

Initial distribution for a straight diamond crystal. The black full curve represents the analytical expression of Eq. (15) of Backe et al. [3], the red dots are simulation calculation with 10,000 trials

7 Results

7.1 Results for channeling at (110) planes of diamond and discussion

Now all ingredients are available to investigate the de-channeling process. Simulations have been performed for a target thickness of 200 \(\mu \)m for three consecutive categories, primary channeling mode, de-channeling mode, and secondary channeling mode. The channeling history is terminated either if the electron leaves the potential boundaries in the secondary channeling mode, if it will not be de-channeled after a distance of 100 \(\mu \)m, or if it reaches the maximum transverse energy of 10 \(u_0\) = 240 eV. Typical examples of trajectories are shown in Appendix C.

Fig. 6
figure 6

Normalized channel number distribution. Channel 0 is the primary channel in which the particle was originally captured, neighbouring ones in which the particle re-channelled are counted from there

Figure 6 shows the channel number distribution for re-channeling. It is interesting to notice that the largest probability for a re-channeling occurs at the direct neighboring channels. This is expected since de- and re-channeling are both governed by the same scattering distributions shown in Fig. 3 which are largest at small energy transfers. The larger the transverse energy becomes, the smaller the probability for re-channeling will be. The electron may re-channel to channel 0 if it returns after de-channeling, i.e. crossing the boundary \(x = \pm d_p/2\), again to channel 0. The integrated probability for re-channeling up to the second neighboring channels amounts to 33%.

Fig. 7
figure 7

a Primary channel occupation and b approximation of the de-channeling rate as function of the penetration depth z, with details shown in the inset. The channel occupation is normalized to unity for \(N_0\) = 10,000 events. The horizontal error band indicates the mean value \(\overline{\lambda (z)}\) = (0.0624 ± 0.0011)/\(\mu \)m as obtained for \(14 \le z/\mu m \le 100\)

Fig. 8
figure 8

Same as Fig. 7 for re-channeled particles. The horizontal error band indicates the mean value \(\overline{\lambda (z)}\) = (0.0641 ± 0.0014)/\(\mu \)m as obtained for \(14 \le z/\mu m \le 100\)

Figure 7 shows in panel (b) an approximation of the instantaneous transition rate

$$\begin{aligned} \lambda _{de}(z) = -\frac{f_{de}'(z)}{f_{de}(z)} = \lim _{\begin{array}{c} \varDelta z\rightarrow 0 \\ N \rightarrow \infty \end{array}} \frac{\varDelta {(N}(z)/N_0)/\varDelta z }{N(z)/N_0} \end{aligned}$$
(36)

as derived from the simulated channel occupation numbers \(\varDelta {(N}(z)/N_0)/\varDelta z\) and the channel occupation \(N(z)/N_0\) depicted in panel (a). For z less than about 14 \(\mu \)m the instantaneous transition rate is larger than the constant mean value for \(z>\) 14 \(\mu \)m, see panel (b). This is a consequence of the fact that the captured electrons need a certain relaxation depth to achieve statistical equilibrium. The structure of the rate for \(z\le 5 \mu \)m is shown in the inset. The constant mean value at statistical equilibrium is \(\overline{\lambda _{de,1}^{plane}(z)}\) = (0.0624 ± 0.0011)/\(\mu \)m, as calculated for the interval \(14 \le z/\mu m \le 100\). A fraction of 18% of particles de-channel in the mean within the relaxation depth.

A de-channeled particle may re-channel. Here such an process is defined by the requirement that the particle must be reflected from the potential wall at least once. Results for these secondary channeling process are shown in Fig. 8. The rate exhibits similar features like that shown for the primary process. The equilibrium depth turns out to be \(\overline{\lambda _{de,2}^{plane}(z)}\) = (0.0641 ± 0.0014)/\(\mu \)m for \(14 \le z/\mu m \le 100\) which is in accord with the primary de-channeling rate. However, a much larger fraction of 50% is lost within the relaxation depth \(z<14\mu \)m.

Alternatively, an effective de-channeling length \(\varLambda _{de}\) may be defined via the integral \(\int _{z_0}^{\varLambda _{de}} \lambda _{de}(z)~dz\) = 1 with \(z_0\) the entrance position into the channel. One obtains from Figs. 7b and 8b \(\varLambda _{de} \approx \) 13 and \(\approx \) 8 \(\mu \)m for the first and second de-channeling length, respectively. These significant different values may be explained by the fact that re-channeling happens mainly in the middle of the channel where the overlap with the atomic density is largest. Therefore, a shorter equilibration time is needed in comparison with the primary channel occupation probability which is distributed over the whole channel width.

7.2 Channeling in bent crystals

For production of undulator-like radiation the crystal must be bent. In the following, the diamond crystal is assumed to have dimensions like the silicon crystal described in Ref. [25], i.e., thickness of 30.5 \(\mu \)m and bending radius of 33.5 mm. Calculations were performed as described in section 7.1. A primary de-channeling rate \(\overline{\lambda _{de,1}^{bent}(z)}\) = (0.0668 ± 0.0010)/\(\mu \)m has been obtained for the interval \(14 \le z/\mu m \le 100\) for 75 % of particles imping the crystal. For re-channeled particles the de-channeling rate turned out to be \(\overline{\lambda _{de,2}^{bent}(z)}\) = (0.0685 ± 0.0025)/\(\mu \)m, for the same interval with an efficiency of only 14 % of particles imping the crystal. While the primary de-channeling rate for the bent crystal increases by only 7 % in comparison with the plane one, the efficiency decreases dramatically.

In addition, the beam profile at the exit of the crystal has been calculated. The result is shown in Fig. 9. The profiles exhibit striking similarities with the results shown by Mazzolari at al. [25, Fig. 3]. In particular the deflection peaks due to channeling in (a) and (b), as well as the volume reflection shift in opposite direction for the tilted crystal in (a) can clearly be recognized.

Fig. 9
figure 9

Beam profiles for a 30.5 \(\mu \)m thick diamond crystal with a bending radius of 33.5 mm for a an entrance angle \(\psi _0\) = -0.45 mrad, and b \(\psi _0\) = 0 mrad

7.3 Comparison with results from the MBN explorer package

Simulation calculations have been performed on the basis of the rather sophisticated MBN Explorer package as described in a number of publications of the group around A.V. Korol, and A.V. Solov’yov et al. In the paper of A.V. Pavlov et al. [7] simulation calculations were done also for channeling of 855 MeV electrons in diamond single crystals. Unfortunately, it is not possible to compare the results directly with results obtained in this paper. Therefore, Fig. 5c (red curve) of Ref. [7] was digitized. The result is depicted in Fig. 10a. In Fig. 10b an approximation of \(\lambda _{p}(z)\) according to Eq. (36) is shown. Particularly striking are the zero-valued de-channeling transition rates for small penetration depths. This behaviour reflects the acceptance criterion, meaning that a particle is accepted to be in the channeling mode only if it changes sign of its transverse velocity twice. The mean equilibrium de-channeling length obtained in the interval between \(10\le z/\mu \text {m} \le 20\) amounts to \(\lambda _{p}\) = 0.083/\(\mu \)m is about (33±16)% larger than that one obtained in this work. Whether this is a significant difference, or not, is in view of the rather large uncertainty currently not clear.

Fig. 10
figure 10

a Digitized fraction \(N_p/N\) of 855 MeV electrons in a straight (110) diamond crystal channel as taken from Pavlov et al. [7, Fig. 5(c), red curve], and b differential quotient normalized to the fraction of electrons in the (110) channel. The estimated equilibrium de-channeling rate in the shown interval \(\overline{\lambda (z)}\) = 0.083/\(\mu \)m has an error of about ±0.01/\(\mu \)m, mainly due to digitization

8 Impact on the Fokker–Planck equation

The solution of the Fokker–Planck equation has been described in our previous work [3] for planar channeling experiments in straight crystals with 855 MeV electrons from the Mainz Microtron MAMI. The drift coefficient, and in turn also the diffusion coefficient, were calculated with the Kitagava–Ohtsuki approximation [6] which in [3] was represented in the following form

$$\begin{aligned} D_1 (E_\bot )= & {} \frac{E_s^2}{2 pv X_0}\times \int \limits _{x_{min}}^{x_{max}}\frac{\textrm{d}P}{\textrm{d}x}(E_{\bot },x) \nonumber \\{} & {} \times \frac{d_p}{\sqrt{2\pi }u_1} \exp (-x^2/2u_1^2)~\textrm{d}x \end{aligned}$$
(37)

with

$$\begin{aligned}{} & {} \frac{\textrm{d}P}{\textrm{d}x}(E_{\perp },x)=\frac{2}{T(E_{\perp })\cdot v}\sqrt{\frac{\gamma \cdot m_e c^2}{{2\left( E_{\perp }-u(x)\right) }}}~~, \end{aligned}$$
(38)
$$\begin{aligned}{} & {} T(E_{\perp })\cdot v = 2\int \limits _{x_{min}}^{x_{max}} \sqrt{\frac{\gamma \cdot m_e c^2}{2\big (E_{\perp }-u(x)\big )}}~~\textrm{d}x. \end{aligned}$$
(39)

The integration limits \(x_{min}\) and \(x_{max}\) are roots of the equation \(E_\bot =u(x)\). The quantity \(X_0\) = 0.1213 \(\cdot 10^6~\mu \)m is the radiation length [26]. The results of the solution of the Fokker–Planck equation depend strongly on the choice of the parameter \(E_s\) since it enters quadratically. In our previous work \(E_s\) = 10.6 MeV was chosen [27], resulting in a de-channeling length of 41.04 \(\mu \)m. This value deviates by a large factor of about 2.6 from the simulation calculation result of 1/\(\overline{\lambda (z)}\) = 16.0 \(\mu \)m quoted in the previous Sect. 7.1. In the following the diffusion coefficients entering the Fokker–Planck equation will be reanalyzed for both, (i) scattering at atoms only, and (ii) including scattering at electrons with the formalism described above.

Fig. 11
figure 11

a Drift coefficients \(D_1\), and b diffusion coefficient \(D_2\) as function of the normalized transverse energy \(E_\perp /U_0\) for a straight diamond crystal with \(U_0 = u_0\). Red curve: contributions from atoms, blue from electrons, and black is the total

Fig. 12
figure 12

Transverse probability distributions. Red full curves depict calculations with Eq. (38), blue dots simulations. a \(\overline{E_\perp }\) = 29.8 eV well above the barrier, b 24.2 eV at the barrier, c 19.7 eV below the barrier, and d 4.7 eV well below the barrier. The shadowed region indicates the nuclear density distribution

(i) For scattering at atoms drift and diffusion coefficients may be written as

$$\begin{aligned}{} & {} D_1^{(at)} (E_\bot ) = \frac{pv}{2}<\theta _x^2>_{at}\frac{8}{a^3}~\cdot \sigma _{tot}^{(at)} \nonumber \\{} & {} \times \int \limits _{x_{min}}^{x_{max}} \frac{\textrm{d}P}{\textrm{d}x}(E_{\bot },x) \frac{d_p}{\sqrt{2\pi }u_1} \exp (-x^2/2u_1^2) \textrm{d}x \end{aligned}$$
(40)

and

$$\begin{aligned} D_2^{(at)} (E_\bot )= & {} \frac{pv}{2}<\theta _x^2>_{at}\frac{8}{a^3}~\cdot \sigma _{tot}^{(at)}\times \int \limits _{x_{min}}^{x_{max}} 2 (E_\bot -u(x)) \nonumber \\{} & {} \times \frac{\textrm{d}P}{\textrm{d}x}(E_{\bot },x) \frac{d_p}{\sqrt{2\pi }u_1} \exp (-x^2/2u_1^2) \textrm{d}x. \end{aligned}$$
(41)

In Eqs. (40) and (41) three terms can be distinguished which have the following meanings: The pre-factor \((pv/2)<\theta _x^2>_{at}\) is the mean transverse energy gain for a single elementary scattering process, the pre-factor \((8/a^3)~\cdot \sigma _{tot}^{(at)}\) represents the mean number of collisions per unit path length, and the integrals represent averages for one oscillation period accounting for the distribution of the scattering centers across the channel. Defining in analogy to Eq. (37) an \(E_{s,m}^2 = (pv)^2 X_0\) \(<\theta _x^2>_{at} \cdot (8/a^3) \cdot \sigma _{tot}^{(at)}\), an microscopically defined scattering parameter \(E_{s,m}\) = 12.13 MeV follows. With this value the de-channeling length as derived from the solution of the Fokker–Planck equation reduces from \(L_{de}^{FP}\) = 41.04 \(\mu \)m for \(E_s\) = 10.6 MeV to the significantly lower value 31.3 \(\mu \)m.

(ii) For scattering at electrons the diffusion coefficients are defined in analogy to the scattering at atoms as

$$\begin{aligned} D_1^{(el)} (E_\bot )= & {} \frac{pv}{2} <\theta _x^2>_{el}\nonumber \\{} & {} \times \int \limits _{x_{min}}^{x_{max}} \frac{\textrm{d}P}{\textrm{d}x}(E_{\bot },x)~n_{el}(x)~\sigma _{tot}^{(el)} \textrm{d}x~~~~~ \end{aligned}$$
(42)

and

$$\begin{aligned} D_2^{(el)} (E_\bot )= & {} \frac{pv}{2} <\theta _x^2>_{el}\int \limits _{x_{min}}^{x_{max}} 2 (E_\bot -U(x)) \nonumber \\{} & {} \times \frac{\textrm{d}P}{\textrm{d}x}(E_{\bot },x)~n_{el}(x)~\sigma _{tot}^{(el)} \textrm{d}x. \end{aligned}$$
(43)

The Fokker–Planck equation has been solved for the sum of the atomic and electronic contributions

$$\begin{aligned}{} & {} D_1(E_\bot ) = D_1^{(at)}(E_\bot )+D_1^{(el)}(E_\bot ), \end{aligned}$$
(44)
$$\begin{aligned}{} & {} D_2(E_\bot ) = D_2^{(at)}(E_\bot )+D_2^{(el)}(E_\bot ). \end{aligned}$$
(45)

Drift and diffusion coefficients are shown in Fig. 11. Although the contribution of the electrons seems to be small, scattering at electrons has a significant impact, i.e., the de-channeling length reduces to 21.8 \(\mu \)m. This value is still 36 % larger than that one obtained by simulation calculations.

In order to examine possible reasons for the remaining difference between the simulation calculation and the solution of the Fokker–Planck equation, in Fig. 12 simulation results for the probability distribution are compared with calculations of \(\textrm{d}P/\textrm{d}x(E_\perp ,x)\) with the aid of Eq. (38). Both probability distributions look very similar. Only a rather small enhancement of the simulation results at the overlap region with the atomic density distribution may be discernable in panel Fig. 12b for \(\overline{E_\perp }\) = 24.2 eV. Tikhomirov [5, ch. 3.4] discussed that the probability distribution Eq. (38) cannot be applied because of large transverse energy variations at domains where the atomic density is large. This statement can more or less be ruled out for diamond, however, such fluctuations may well appear for crystals composed of atoms with large atomic numbers Z.

9 Discussion and conclusions

The instantaneous transition rate concept has been applied for channeling of particles in straight and bent crystals. It turns out that the de-channeling length of electrons in single crystals cannot be described by a single number since the transition rate approaches a constant value only after a certain relaxation length. It is in the order of about 15 \(\mu \)m for our example of 855 MeV electrons channeling in (110) planes of diamond. A reasonable approach to characterize a de-channeling length would be the assignment of two numbers, the equilibrium de-channeling length and the fraction of particles for which it holds. As demonstrated above, these fractions drastically differ for the primary and secondary channeling processes.

The main deficiency of the Fokker–Planck equation may be found in the elimination of the transverse x coordinate by averaging out any details associated with the dynamics and the distribution of the scattering centers across the transverse channeling coordinate. This fact reveals itself in the instantaneous drift and diffusion coefficients for each step (\(\varDelta x_k,\varDelta z_k \)), see Eqs. (34) and (35). The mentioned deficiency can probably not be cured by a parameter adaption. The Fokker–Planck equation can, therefore, only be used to get a qualitative idea of the channeling process.