Simulation code for modeling of coherent effects of radiation generation in oriented crystals

We present the CRYSTALRAD simulation code, combining all the features of the CRYSTAL simulation code for simulations of charged particles trajectories in a bent crystal and the RADCHARM ++ code for calculation of the radiation spectrum. The CRYSTALRAD code is based on Monte Carlo simulations of trajectories in the planar and axial electric field either in a straight, bent, or periodically bent crystal taking into account multiple and single Coulomb scattering on nuclei and electrons, nuclear scattering and ionization energy losses. The trajectories simulated are used for calculation of radiation spectra by the Baier-Katkov method. We compare our simulations with experimental data taken at MAMI (MAinzer MIkrotron) as well as give an example for a possible future study with sub-GeV electrons interacting with Si bent crystals.


I. INTRODUCTION
Since the 1950s, it has been known that the lattice structure can strongly influence the electromagnetic processes in oriented crystals and that the alignment of a beam of electrons/positrons with respect to the axis or planes of a crystal leads to an increase of the probability of bremsstrahlung emission.
Coherent orientational effects in a crystal can be exploited for various applications in accelerator physics as well as for the development of novel x-and gamma-ray crystal-based sources. In particular, coherent bremsstrahlung (CB) [1] is currently exploited in several laboratories worldwide, in which electron beams of energy of the order of GeVare available, such as MAMI (Mainz, Germany) [2], JLAB (Newport News, USA) [3], ELSA (Germany), and Max Lab (Lund, Sweden), to generate intense monochromatic and linearly polarized gamma beams for experiments in photonuclear and hadronic physics researches. In addition, at certain orientation of a scintillation crystal it is possible to accelerate the electromagnetic shower development, which is very attractive for forward e.m. (electromagnetic) calorimeters technology [4][5][6]. Furthermore, bent crystals have already shown themselves as cheap and efficient tool for crystal-based manipulation, for application in hadron beam collimation and extraction at U70 at IHEP-Protvino, Tevatron at Fermilab and SPS and LHC at CERN [7][8][9][10][11][12].
Coherent interaction of charged particles with crystals occur when the particle trajectory is oriented with a small angle with respect to a main crystal axis or plane. In such a case, correlated collisions with atoms in the same row/plane occur and therefore the charged particle interacts with crystalline planes or atomic strings (axes) as a whole. This modifies not only the motion of charged particles inside crystals, but also the process of bremsstrahlung radiation emission. First the increase of the bremsstrahlung in oriented crystals was predicted to appear due to an effect called coherent bremsstrahlung (CB). CB was then discovered in Frascati in 1960 and consists of the enhancement of bremsstrahlung when the momentum transferred by electrons to the crystal matches a reciprocal lattice vector, in analogy with the Bragg/Laue diffraction [1]. When the charged particle velocity is nearly parallel to either a crystal axis or a plane, the particle trajectory would be forced in an oscillatory motion within the planar/axial potential well, i.e., the channeling phenomenon occurs, leading to a specific e.m. radiation emission called channeling radiation (CR) [13]. Channeling may occur if the incidence angle of the charged particle with respect to the crystal planes/axes is lower than the critical value introduced by Lindhard θ c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2V 0 =εÞ p , where ε is the particle energy and V 0 the average potential energy amplitude [14].
Channeling in bent [15][16][17][18][19][20] and periodically bent crystals, the latter case corresponding to a crystalline undulator [21][22][23], can be applied for the generation of x-ray and gamma radiation. Monte Carlo simulations are needed to predict the characteristic of the emitted radiation and need to combine the simulation of electron/positron motion in crystalline materials as well as the generation of radiation, which is not trivial. Indeed, the particle motion includes different types of scattering in an ordered crystalline media, the cross section of which strongly depends on a particle transverse position, while the radiation generation at high energies cannot be described in terms of classical electrodynamics because of low probability and high energies of photons produced.
There are different approaches of trajectory calculation in a crystal, which are reasonable for application at different energies. At low energies, i.e., few MeV for electrons, only quantum mechanics allows one to treat the channeling effect correctly [24][25][26][27], while starting from the energies of ∼10 2 MeV a classical description works well. In the classical approach, one may utilize the scheme of binary collisions which assumes that the motion of a projectile at all times is influenced by the force due to the nearest atom or, with the sufficient computer power, one may account for the interaction with larger number of the crystal atoms [28][29][30][31][32]. The limitations of the classical binary collision model are discussed in the work [33].
Another possible approach in the classical framework is to simulate the particle trajectory in the average continuous planar or axial potential, proposed by Lindhard [14]. In this model, already realized in a number of computer codes [34][35][36][37][38][39], one simulation step usually contains hundreds of atoms, considerably reducing the calculation time. The main difference between these computer codes is the realization of atomic potential as well as the Coulomb scattering model, strongly affecting charge particle trajectories. The main approach of radiation spectra simulations in the energy range, starting from hundreds of MeV is the Baier-Katkov (BK) quasiclassical operator method [28], in which the probability of radiation emission is computed through integration on the classical particle trajectory, while at the same time it takes into account the quantum recoil of the primary particle in the photon emission. Thereby the BK method is suitable with the models in which the particle motion is described by classical trajectories, in particular for the Lindhard averaged potential one. In the literature, different algorithms are present for the integration of the radiation emission probability in the framework of the Baier-Katkov method realized [16,32,[40][41][42], with differences in the integration methods and, consequently, also in calculation speed.
In this paper, we present the CRYSTALRAD simulation code providing fast Monte Carlo simulations of both charged particle dynamics and radiation emission in straight, bent and periodically bent crystal of any material and crystal lattice type with well verified experimentally models of scattering and radiation [16][17][18][19]34,40,[43][44][45]. Apart from this, the CRYSTALRAD code includes a wide number of features, namely peculiarities of crystal geometry, possibility of simulations for a set of initial parameters to solve the optimization problem and, furthermore, MPI parallelization provides the opportunity to carry out complex calculations at supercomputers with a huge reduction of needed simulation time.

II. PHYSICAL MODEL FOR TRAJECTORIES CALCULATION
A. Trajectories calculation in the averaged continuous atomic potential A trajectory of charged particle, moving at relatively small angle with respect to the crystal planes or atomic strings can be calculated in the approximation of the averaged atomic potential, as first proposed by Jens Lindhard [14]. In other words, one needs an averaged interplanar or interaxial potential UðxÞ to integrate the trajectory equations, that can be written as where p is the particle momentum, v its velocity, x and y transverse as well as z longitudinal coordinates, U 0 x and U 0 y components of transverse electric field, R x and R y horizontal and vertical bending radius of a bent crystal (in the case of a straight crystal 1=R x ¼ 1=R y ¼ 0). In the approximation of the interplanar potential one may choose the direction of y-coordinate parallel to the crystal plane, which causes U 0 y ¼ 0, and consequently disappearance of the second equation in Eqs. (1).
There are a lot of single atom potential models, most of them are defined through atomic form factors measured with x-ray diffraction [46]. The examples of interplanar and interaxial potential, obtained for the Doyle-Terner approximation [47] with the coefficients calculated in [48] in comparison with harmonic and Molière approximations [49] are presented in Fig. 1.

B. Scattering model
The solution of Eqs. (1) gives only coherent part of particle-crystal interaction, but does not include the incoherent one, i.e., multiple and single scattering on separate atoms, nuclei and electrons, for instance Coulomb scattering. A particle can be scattered by a screened atomic potential as a whole or by a single electron, causing atom ionization.
Coulomb scattering on atom can be treated by a well fitted cross section of screened atom [50][51][52], based on Yukawa potential: where ϑ ¼ , ϑ x and ϑ y are scattering angles in the horizontal and vertical planes, respectively, acteristic scattering angle [50,51], z p particle charge, Z atomic number of atom, e electron charge, a TF ¼ ð9π 2 =128ZÞ 1=3 a 0 Tomas-Fermi screening radius, a 0 Bohr radius, dΩ an element of a solid angle. By using this cross section and taking into account dΩ ¼ ϑdϑdφ, where φ is the polar angle, one can find the r.m.s. angle of multiple scattering in a material of thickness Δz. However, one must take into account that in a crystal the scattering cross section is divided onto coherent and incoherent part [53], where the incoherent part can be expressed by using the Debye-Waller factor expð−p 2 ϑ 2 u 2 1 Þ. Therefore the r.m.s. multiple scattering angle can be written as: where n N is the average nuclear density along the trajectory element Δz. The angle ϑ 2 is the angle, splitting multiple and single scattering [53]. Indeed, according to Molière [50,51], Coulomb scattering involves both multiple scattering, described by a Gaussian distribution, and single scattering characterizing the long non-Gaussian tails of the distribution of scattering angles. ϑ 2 is not strictly defined, since at its zero value only single scattering is actually left. However, the scattering on small angles can be taken into account as a multiple scattering considerably reducing the number of necessary computations. One can distinguish between single and multiple scattering by calculating ϑ 2 using the probability of a single scattering. For convenience of modeling, the value of probability equal to 0.1 was used. The length at which the next Coulomb occurs single scattering can be written as: Single scattering suppression is taken into account with cutoff by specially generated random number, compared to the value of the Debye-Waller factor [53].
Another type of single scattering event is scattering on electrons. Its cross section is represented by the Rutherford cross section [54]: However, the main parameter, defining scattering on electrons is the kinetic energy, transferred to electron T e ¼ p 2 ϑ 2 =2m e , where m e is the electron mass, which can be introduced in (5) during integration. Consequently, the length of the next single Coulomb scattering on electron can be calculated similarly to (4): where n e is the electron density, T max ¼ 2m e ðγβÞ 2 = ð1 þ 2m e γ=m þ ðm e =mÞ 2 Þ the maximal energy that can be lost for ionization [55], while T min is the ionization potential that can be found as the Fermi energy at the point of scattering, β ¼ v=c.
Nuclear scattering is important only for channeling hadrons (protons, ions). The length of nuclear scattering events L el , L diff , L inel (elastic, diffractive, and inelastic, respectively) is almost independent of the energy and can be usually found in the literature [55] for any type of material. One must only recalculate this length, taking into account the nuclear density along the trajectory: The r.m.s. scattering angles can be written as: where B is a slope determined for each type of scattering. For example, for protons in the case of elastic scattering [56], it is written as: where A is the atomic number of a nucleus in a crystal. All the types of scattering depend on nuclear and electron density distributions. The examples of these distributions are shown in Fig. 2. One can notice that the electron density for germanium is more than twice higher, resulting in higher probability of single scattering on electrons than for Si.

III. ALGORITHM FOR THE DIRECT INTEGRATION OF THE BAIER-KATKOV FORMULA TO COMPUTE THE RADIATION EMISSION
The BK method is based on classical trajectories of ultrarelativistic particles. Therefore, a QED problem is reduced to the solution of classical equations of motion of a charged particle in crystalline field by integration along the classical trajectory. Moreover, this method takes into account the quantum recoil that is not negligible at ultrarelativistic energies and thereby is applicable in the whole energy range for photon energies ω, except the extreme case when photon energy approaches to the particle energy, i.e., ω ≅ ε (we are working a unit system for which The BK quasiclassical formula [28] can be written as: where N 21 is the radiation polarization matrix, α is the fine structure constant, k ¼ ðω; kÞ is the 4-vector of the photon momentum including radiated energy ω and 3-momentum k, k 0 ¼ εk=ε 0 , where ε and ε 0 ¼ ε − ω are the particle energy before and after the photon emission, respectively, x 1;2 are the particle coordinate 4-vector. In the approximation of small angles, Eq. (10) can be written as [16,40]: γ is the Lorentz factor. Thereby, the evaluation of the radiated energy (1) is reduced to the calculation of the integrals jI ⊥ j 2 and jJj 2 . In order to accelerate the numerical integration, one simplifies the integrals through an integration by parts [16,40]: where ϕ ¼ k 0 xðtÞ,t i ¼ ðt i−1 þ t i Þ=2 and N is the total number of steps. A similar procedure can be implemented for the integral I. The result of integration leads to the radiation spectral intensity dE=dω, namely the value that is usually measured experimentally.
The Baier-Katkov method is specially designed for the cases, when the energy of radiated photon is comparable with the particle energy. For other cases such as transition radiation or Cherenkov radiation, calculation approaches based on classical electrodynamics can be applied. The equations of classical electrodynamics are also applicable to the cases of radiation in photonic crystals [57,58] and 2D materials [59].

A. General description
The CRYSTALRAD simulation code is an unification of the CRYSTAL [34] simulation code and the RADCHARM++ [40] routine. The CRYSTAL simulation code is designed for trajectory calculations taking into account various coherent effects of the interaction of charged particles (protons, ions, electrons, muons, pi-mesons, and their antiparticles) with straight or curved single crystals and different types of scattering as described in Sec. II. The RADCHARM++ code contains the routine to calculate single-photon radiation spectra through the integration of the BK formula with the algorithm presented in Sec. III, by using of charged particles trajectories, calculated by the CRYSTAL code.
The program contains one-dimensional and twodimensional models that allow one modeling of classical trajectories of relativistic and ultrarelativistic charged particles in the field of atomic planes or strings, respectively. The energy range accessible by the CRYSTALRAD code is from 100 MeV for electrons/positrons as well as 1 MeV and higher for protons and ions. Indeed, in case of lower energy for electrons or positrons quantum effects can contribute significantly and the classical trajectory approximation fails to work. The restriction on the angles of particle incidence with respect to the crystalline axes or planes is no more than several degrees. Otherwise, the longitudinal step becomes substantially smaller than the interatomic distance, and both one-dimensional and twodimensional models are inapplicable. For calculation of radiation spectra the crystal length is limited by the requirement of single-photon emission, given by the total probability of photon emission to be less than few tenths.
The algorithm for simulation of motion of particles in a crystal is shown in Fig. 3. The trajectory of a charged particle is calculated by numerical solution of the equations of motion (1) by the fourth-order Runge-Kutta method "3/8" rule [60,61], which gives the dependence xðzÞ, yðzÞ, θ x ðzÞ, and θ y ðzÞ, where θ x ðzÞ ¼ dx=dz and θ y ðzÞ ¼ dy=dz are the particle angles in a horizontal and vertical plane, respectively.
Coulomb scattering is modeled taking into account the suppression of incoherent scattering by the model [53] (2)-(3), as described above. The simulated scattering angles, denoted by ϑ x , ϑ y in Fig. 3, are added to the angles of the particle, obtained by solving the equation of motion. In addition to this, nuclear elastic, diffractive, and inelastic scattering are simulated by the model [56] according to Eqs. (7)- (9).
Geometry of a crystal is also included in the program, namely crystal bending, the possibility of entry/exit through its lateral surface, the influence of the miscut angle [62,63], as well as a crystalline cut, which significantly increases channeling efficiency [64,65]. A special procedure allows one to vary one or more parameters by specifying a table of values for them and performing calculation for each set of parameters. Output files include the values of horizontal and vertical coordinates and particle angles at the entrance and exit of the crystal, as well as the efficiency of channeling, the number of inelastic scattering events, etc.
For the calculation of trajectories and scattering, the following functions are used: interplanar potential, interplanar electrostatic field, nuclear and electron densities, minimum ionization energy of an atom. All these functions are stored as interpolation coefficients of cubic splines [60,61] in the input file, which allows one-dimensional and two-dimensional models to be specified for any monocrystal of any orientation without modifying the source code. The interplanar potential and its derivative functions were preliminarily calculated by using the Doyle-Turner potential [48].

B. Implementation of MPI parallelization
One of the main advantage of the CRYSTALRAD simulation code is the application of the message passing interface (MPI) allowing one to carry out parallel calculations on clusters and supercomputers. Since the trajectory of each particle is simulated independently, MPI parallelization does not require a lot of communications. With the help of the procedures MPI_COMM_SIZE and MPI_ COMM_RANK, each process receives equal parts of the particles for simulations. The uniqueness of each particle trajectory is ensured by the individual initialization of the random number generator for each process, in which the process number is the argument. Since MPI is designed for distributed memory systems, all the input data (interpolation coefficients of the main functions, etc.) are loaded into the RAM separately for each process. In particular the first process reads all the input data and then broadcasts these data to the other processes by using the MPI_BCAST procedure.
The main limitations of linear scaling of performance is the input procedure. However, since it is carried out only once at the start of the program run, it does not contribute considerably. Moreover, since simulations do not practically require interactions between different processes, the productivity indeed increases almost linearly with their number. Data exchange is performed only at the end of the calculations by using the MPI_Reduce procedure, which calculates the channeling efficiency, the fraction of inelastic losses, and other data that are written to the output file.
Since it concerns only a few numbers, it practically has no effect on performance.

C. CRYSTAL and RADCHARM++ unification
The code RADCHARM++ was included into the CRYSTAL code as a subroutine. The input parameters of RADCHARM++ are arrays containing charged particle angles θ x and θ y as well as scattering angles ϑ x and ϑ y for each trajectory step. In addition the spectra parameters, namely spectrum energy range and step are also uploaded to RADCHARM++. The output result is the spectrum for each trajectory, being transferred as an array to the CRYSTAL code and then averaged over all the trajectories.
In order to verify the validity of the CRYSTALRAD code, the radiation spectra in the experiment [19] at Mainz Mikrotron MAMI have been simulated for different alignment of a silicon bent crystal. The experimental setup at MAMI is well described in papers [19,43,66]. In the experiment, a 855 MeV electron beam with low divergence, i.e., considerably less than the critical channeling angle, was directed through a bent crystal, aligned along its (111) bent crystal planes. After interaction with the bent crystal, charged particles were deflected by a magnetic field and then registered by Si microstrip detectors [19,43,66]. Photons produced inside the crystal and moving straightforward were detected by a NaI scintillator detector coupled to standard PMTs. The crystal was 30.5 μm thick with a bending radius of 33.5 mm. The angular divergence of 855 MeV electron beam was 30 μrad. The simulation results in comparison with the experimental data are presented in Fig. 4, confirming the correctness of our simulations. In both simulations and experiment, the crystal was aligned at 0 and 493 μrad with respect to the beam direction, corresponding to channeling and volume reflection orientations, respectively. Volume reflection is the deflection of overbarrier particles occurring in a bent crystal, when particle trajectory becomes parallel to crystal planes in the crystal volume due to the bending itself [39].
Furthermore, we simulated the expected radiation spectra for experimental cases presented in [45,67]. The scope of the experiment [45,67] was the measurement of channeling characteristics, i.e., channeling efficiency and dechanneling length, in dependence on the radius of bent crystals. Since channeling radiation strongly depends on these characteristics, a study of channeling radiation spectra vs crystal radius could be relevant for the development of novel x-ray and gamma radiation sources. This dependence has been obtained by using CRYSTALRAD simulation code for the silicon sample. The experimental setup was the same as in the experiment [66] with the replacement of the Si microstrip detector with a LYSO screen. In details, a 855 MeV electron beam impinged on a 15 μm thick Si crystal bent along its (111) planes. The angular divergence of the beam was 21.4 μrad, being one order on magnitude less than the critical channeling angle. The simulated spectra in dependence on the radius of curvature of the crystal, divided by the critical channeling radius R cr ¼ E 0 =pv, where E 0 is the maximal interplanar electric field, are presented in Fig. 5. The critical radius R cr represents the minimal radius of curvature to still have channeling.
In Fig. 5, one may notice the CR peaks at the energy of the order of 2 MeV (as in Fig. 4), while the radiation intensity increases with the radius of curvature, R. Such behavior would be interesting to investigate in future experiments.

V. CONCLUSIONS
The CRYSTALRAD simulation code has been introduced. It is designed for simulation of coherent effects of charged particles interaction with crystals, accompanied by radiation generation. The code is based on classical trajectories simulations and calculation of radiation spectra by the Baier-Katkov quasiclassical method. The code includes a wide number of features, allowing one to carry out simulations in crystals of any material and any crystal lattice type and shape, namely can be used to simulate straight, bent and periodically bent crystals. Moreover, the option of variation of initial parameters during one simulation run is included in the code, to reveal the dependence of the system on initial parameters as well as to solve the optimization problem. In addition, MPI parallelization allows one to apply the simulation code on supercomputers. The limitations of the code are connected with the minimal energy of particles to be 100 MeV for electrons/ positrons and with the constraint of spectrum calculation being the single photon approximation, limiting the crystal thickness by the radiation probability to be lower than few tenths.
The CRYSTALRAD simulation code is the unification of CRYSTAL and RADCHARM++ codes, already validated separately in past experiments, with improvement due to MPI parallelization. An example of validation of CRYSTALRAD simulated radiation spectra has been provided for the experiment at Mainz Mikrotron MAMI where 855 MeV electrons interacted with a 30 microns bent Si crystal. Moreover, simulations of a possible experiment to study channeling radiation vs. curvature radius at the Mainz Microtron MAMI have been carried out for a silicon bent crystal of 15 microns thickness. The radiation intensity increases with the radius of curvature. This would be interesting to investigate in view of application for a novel type of gamma-ray source based on bent [19] or periodically bent [21][22][23] crystal.
In conclusion, the CRYSTALRAD simulation code can be applied for the simulations of a wide number of applications of crystal for both beam steering and radiation in a wide range of the energies of modern accelerators, from hundreds of MeV up to tens of TeV.