Theoretical and Simulation Studies of Transverse Beam Size Effects on Optical Transition Radiation in Prewave Zone

Optical transition radiation (OTR) has been studied and applied on the beam diagnostics for decades. The potential implication of OTR also includes THz radiation sources. Therefore, the theoretical analysis and simulation tools have become indispensable for the OTR study. Moreover, the OTR theory, for the wave zone (far-field approximation), has been widely used in the literature. However, these theories for the prewave zone have been proposed on the basis of single electron approximation. In this study, we developed a theory with consideration of the electron beam structure based on Kirchhoff’s method for studying the effect of beam transverse size on the angular distribution of the OTR in prewave zone. The proposed formalism involves complicated con-volution integral of functions and the dimensions of integrand are not low. To perform such integral, we developed a Fortran program for quasi-monte Carlo method, which is robust and suitable for high dimensional integration. The disadvantage of this method is that a large amount of samplings may need to be employed to achieve good convergence. Therefore, in order to get the radiation angular profile, we need to perform such integration for different observation angles that might be computationally intensive. Therefore, we parallelize the program with the Message Passing Interface (MPI) programming concepts. To verify the theoretical calculations, few two-dimensional FDTD simulations were carried out that show the validity of the proposed model. The proposed theory and numerical tool would be used to predict radiation properties of the NSRRC THz coherent transition radiation (CTR) in the prewave zone.


Introduction
When an electron traverses across an interface of two different media, radiation is emitted.
is phenomenon is called transition radiation (TR) and was theoretically predicted by Ginzburg and Frank in 1946. eir approach to tackle this problem is based on solving the electromagnetic wave equations directly with consideration of boundary conditions at the interface [1][2][3] and provided solutions in wave zone (far-eld). It is often called optical transition radiation (OTR). If the wavelength is not shorter than that of the visible light, OTR from metallic foils have been used for electron beam diagnostics for decades [4][5][6]. In such applications, metallic foils are considered as perfect conductors and the Ginzburg-Frank formula has been widely used for the theoretical prediction of the radiation properties. Alternative to their approach, as suggested by Fermi [7], the velocity elds of the incident electron are decomposed into Fourier components of waves such that the problem can be modeled as scattering of waves by the interface. Kirschho 's method can be used to solve this scattering problem [8,9]. With this integral formalism, one may nd that the condition for far-eld approximation holds when the observation distance r 0 >> c 2 λ. In this regard, the radiation source is the surface current on conductor which is induced by the electric eld from incident electron; the size, or more specifically the cut-off size, is of the order c λ because the asymptotic behavior of particle field falls off with radial distance d perpendicular to the trajectory that is proportional to exp(− 2πd/βcλ). Due to the finite source size, the prewave zone can be defined as c 2 λ > r 0 >> cλ. e behaviors of the OTR in this region, such as spectral angular distribution, are quite different from those in wave zone [10,11]. e OTR properties in the prewave zone attracts much attention in applications, such as beam profile monitoring and generation of THz radiation in high energy accelerators because the far-field condition may not be satisfied [12,13]. However, far-field approximation is usually used in the experimental designs [14]. On the other hand, as pointed out by Orlandi et al., beam transverse size would affect the behavior of the angular spectrum and may be an important issue for beam diagnostic. In their work, the authors stated that the bunch size effects for normal incident beam of OTR in far-field were studied [15][16][17][18]. For a bunch with Gaussian distribution exp(− (x 2 + y 2 /2σ 2 ⊥ ))exp(− (z 2 /2σ 2 ‖ )), the OTR angular spectrum is proportional to F ⊥ (λ, θ)F ‖ (λ)I e (λ, θ)A, where I e is the OTR angular spectrum of single particle with normal incidence. Moreover, F ⊥ (λ, θ) and F ‖ (λ) are transverse and longitudinal form factors, respectively. We intend to design a theory with consideration of the electron beam structure based on Kirchhoff's method for studying the effect of beam transverse size on the angular distribution of the OTR in prewave zone. e proposed formalism involves complicated convolution integral of functions and the dimensions of integrand are not low. To perform such integral, we developed a Fortran program for quasi-monte Carlo method, which is robust and suitable for high dimensional integration.
In this research, based on Kirchhoff's method, a theoretical formalism for the prewave zone OTR is proposed to study the beam size effects. A two-dimensional numerical simulation code based on finite difference time domain (FDTD) method is used to verify the theoretical predictions. e theory would be used to predict radiation properties of the NSRRC THz coherent transition radiation (CTR) in the prewave zone [11]. We will review Kirchhoff's method, and the pseudo-photons from a bunch of electrons are also introduced and illustrated in detail. After that, a theory for prewave OTR from a bunch of electrons with oblique incidence is derived. Basic concepts about FDTD are introduced, and the simulation model to verify the theory is described. e effects of beam transverse size on prewave zone OTR are discussed, and the simulation verifications of the case with normal incidence are presented as well (using a case study). e major contributions of this research are as follows.
(i) We developed a theory with consideration of the electron beam structure based on Kirchhoff's method for studying the effect of beam transverse size on the angular distribution of OTR in prewave zone.
(ii) We developed a Fortran program for quasi-monte Carlo method, which is robust and suitable for high dimensional integration.
(iii) e developed Fortran program is then used to tackle the complex numerical integration problem.
(iv) We parallelize the program with MPI programming interface. e rest of the paper is organized as follows. In Section 2, we offer an overview of the materials, methods, and theories. Section 3 is about the datasets and evaluation metrics using simulations. Moreover, experimental details are also presented. In Section 4, results are discussed. Moreover, a case study is also discussed. Finally, Section 5 concludes this paper and offers several directions for further research and investigation.

Related Work, Theories, and Models
In this chapter, we will first review Kirchhoff's method for plane conducting target. Based on superposition principle, the velocity fields form a bunch of electrons which are discussed. With velocity fields from a bunch of electrons, a formalism for prewave zone OTR from electron bunch with oblique incidence is derived. e beam size effects in prewave zone and wave zone are also discussed [19].

Kirchhoff's Method for OTR.
In this section, we would review the theoretical treatment proposed by Karlovets [20].
e coordinate system in original work is in left-hand coordinate system and we will follow this convention.

Kirchhoff's Method for Plane
Conductor. Consider an OTR system (Figure 1), in which the perfect conducting plane target is situated at z � 0, and an electron bunch moves with an incident angle α relative to z-axis. To construct the scattering theory for such system, we can start with Green's identity.
where f and g are smooth function, V is a volume in space, S is the boundary enclosing V, and n is the unit direction normal to S. By the method of image, a suitable choice of Green's function for this system is determined through using the following equation: where r 0 ′ is the mirror image of r 0 with respect to plane at z � 0. If we put f � E R i we can get the fields relation for each component of E R . To apply the integral equation to scattering problem, we can make E R (r 0 , ω) � E(r 0 , ω)− E 0 (r 0 , ω) and this yields where E and E 0 are total field and incident field. With the boundary conditions of perfect conducting plane, E x,y (r, ω)| S � 0, we can further conclude 2 Scientific Programming For r 0 >> r, we have and therefore the radiation fields (scattered field) become as illustrated by

Radiation Field in Observer's
Frame. An electron bunch moves with an incident angle α relative to z-axis ( Figure 1). Because OTR can be considered as the particle's fields scattered by the target, it may obey the optical law of reflection. erefore, one should put the center of detector (COD) with angle α relative to z-axis. e position of observer can be described with θ x and θ y relative to COD, and the coordinates transformation between observer's frame and target's frame (k) can be characterized by the two consecutive rotations, as given by and the radiation fields transformation between these two frames is illustrated as Assuming that the longitudinal component in the observer's frame is absent, we get With the fields in observer's frame, the radiation energy spectrum can be calculated with the following equation:

Velocity Field from a Bunch of Electrons.
When investigating the interaction between matter and electrically charged particle, which has a uniform velocity, one may consider this problem as scattering. e moving fields from the charge particle are scattered by the media. is concept was proposed by Fermi [21]. When tackling a scattering problem, one normally considers the incident fields in frequency domain. erefore, we can transform the velocity fields into the frequency domain through the well-known Fourier decomposition method [21]. For the velocity fields from single electron moving along + z e direction, we have When an electron moves in z e direction in space with constant velocity, then there exists a time t 0 , such that this electron would arrive the plane z e � 0 when t � t 0 . If this electron has a deviation in transverse coordinate, say for example (x 0 , y 0 ), the velocity fields in the time domain can be written as given by For a bunch of electrons moving along z e direction with constant velocity v having the distribution NB where N is the number of electrons in this bunch, B ⊥ (x ′ , y ′ ) and B ‖ (t ′ ) are longitudinal and transverse distribution, respectively. It should be noted that the velocity fields from this bunch are, based on superposition of the fields from each electrons, illustrated through the followingequation: where F ‖ � B ‖ (t ′ )exp(iωt ′ )dt ′ is the longitudinal form factor for this bunch. e incident fields in the scattering formalism are represented in the target frame, and we may want to express the fields from electron bunch (the incident source) in such frame. For the OTR system, we consider (refer to Figure 1), the transformation between electron's frame and target's frame is and velocity fields from this electron bunch in the target's frame are

OTR from a Bunch of Electrons.
To formulate the problem of OTR from electron bunch, we begin with equation (6) and use the velocity fields from electron bunch (equation (16)) as the incident field. is yields the radiation fields in the target's frame Transforming to observer's frame by using equation (9), we get 4 Scientific Programming e following formula can be calculated by equations (11) and (16): With the radiation fields, one can calculate the radiation energy with (10).

Numerical Integration.
In the numerical calculation, an electron bunch with the Gaussian distribution is considered, and we calculate the radiation fields in the plane of incidence. In this case and transforming the variables to polar coordinate, (19) becomes ρ 2 � x 2 + y 2 cos 2 α and ρ ′ � x ′ 2 + y ′ 2 ; this yields . To get the angular profile of the prewave OTR from electron bunch, one needs to evaluate equation (19). However, the integrands in equation (19) involve complicated convolution of functions and are not low dimensional. It might be difficult to find the interpolating functions which are easy to evaluate. To tackle this numerical integration, we use Monte Carlo method, which is robust and suitable for high dimensional integrations. e disadvantage of this method is that a large amount of points may, in general, be needed to achieve good convergence. We developed a code with Fortran. In this code, instead of pseudo-random sequence, Sobol sequence is implemented, which is of low discrepancy and can lead to better rate of convergence [22]. To get the angular profile, we need to perform such integration for different angles. is is a time-consuming task, and therefore we parallelize this implementation code with the MPI programming interface and model. is will reduce the time needed to reach particular decisions. e detailed description for the algorithm and source code are given later in the appendices. Scientific Programming 5

Simulations and Experiments
To verify the results from theoretical calculation, simulations based on finite difference time domain method (FDTD) are performed. In this paper, concepts of FDTD method are surveyed, and the simulation model based on FDTD for OTR is described after that.

Finite Difference Time Domain Method.
e FDTD method is a widespread method for computational electromagnetism and was proposed by Yee. We would review some key concepts of FDTD method for two-dimensional system [23]. In isotropic medium, the time evolution of electromagnetic fields can be fully characterized by Maxwell's equations and constitutive relations.
For two-dimensional system with the coordinate (x, y), the fields behavior can be split into transverse electric (TE z ) mode and transverse magnetic (TM z ) mode One may then use central difference scheme to discretize Maxwell's equation and get the following discretization for TE z mode and for TM z mode ( Figure 2).
In most cases, the electromagnetic problems are considered in unbounded regions. It might mean that we have to generate the simulation domain as large as possible to avoid the finite boundary effect. However, this is inefficient or even impossible in practice. One of the most common ways is to create a lossy medium outside the physical domain, and the waves are damped inside this medium. Besides, we also require that the absorbing layer to be reflectionless to the incident waves with all kinds of incident angles. is is the so-called perfectly matched layer. In the FDTD code, the uniaxial perfectly matched layer (UPML) is implemented. e layer is based on the anisotropic medium [24], and the interface between this medium and free space can be reflectionless with suitable choice of the material properties. In this subsection, we will follow and go through briefly the construction of this medium.
Consider a two-dimensional space is divided into region 1 (x < 0) and region 2 (x > 0). Region 1 is isotropic space with the propagating incident wave which propagates toward region 2. e region 2 is full of anisotropic medium which is characterized by e Maxwell equations in this medium are We consider here the case for TE z . To calculate the reflection coefficient, we write down the plane wave in both regions where Γ and T are reflection and transmission coefficient. By using the tangential boundary conditions we can solve the reflection and transmission coefficient 6 Scientific Programming and also get Snell's law Substituting into reflection coefficient, we can see that the interface is reflectionless for any angles of incidence wave.
So far we have discussed the UPML, which is reflectionless for the boundary normal to x direction and attenuates the wave propagating in x direction. We also need to consider the case for y direction. In general, Maxwell's equation in UPML can be expressed as where s is For two-dimension case, we have three kinds of material properties for UPML ( Figure 3

Simulation for OTR with Normal Incidence.
Although the theoretical model developed in previous chapter can be applied for the electron bunch with oblique incidence, we only perform the simulation for the electron bunch with normal incidence. In this case, the spectral angular distribution of OTR is azimuthal symmetry and the simulation  model can be reduced to two dimensions. Our simulation model is based on two-dimensional FDTD. e code we use in this study is written in Fortran with MPI parallelization and is developed by Dr. Toseo Moritaka. e simulation setup is showed in Figure 4.
Once the current pulse moves across the PEC boundary, OTR will be generated and then propagates toward the right hand side. e time signals for H z fields situated at semicircular arc will be recorded. By using FFT, we can transform the time signals for each angle to specific frequency. e absolute square of the fields from different angles in frequency domain is corresponding to the angular distribution of radiation.  8 Scientific Programming normal incidence. We also compare the angular distributions in wave zone with the angular distributions in prewave zone.

Wave Zone vs. Prewave Zone.
We studied the angular distributions from various transverse beam sizes at different observation distances. e outcomes are shown in Figure 5.
We can firstly observe that the angular distributions change significantly with distance and become constant when the distance is large enough, such as r 0 � 5 × 10 3 λ and r 0 � 1 × 10 4 λ. Besides, we can also observe that in wave zone the larger the beam transverse size is, the narrower the angular distribution is. However, this trend can not happen in the prewave zone. Taking angular distributions at r 0 � 1 × 10 2 λ, for example, when the beam transverse size is increasing, the distribution will firstly narrow down and then broaden with the passage of time. is shows that the angular distributions in prewave zone are more sensitive to transverse beam size.

Simulations.
To verify the theory, few two-dimensional FDTD simulations were performed. is should be remembered that in this study: we only verified OTR from electron bunch with normal incidence. In this case, the radiation angular distribution is azimuthal symmetry, and the simulation model can be reduced to two dimensions ( Figure 6).

Simulations vs.
eory. To verify the theory, we performed the simulation for bunch with normal incidence and compared it with the results from theory. e obtained results and findings are shown in Figure 7.
We can see the results from simulation qualitatively match the prediction form theory: in the prewave zone, say r 0 � 1 × 10 2 λ, the angular distribution from larger beam transverse size can be broader than that from smaller beam transverse size. Also, as the distance increases, like Scientific Programming r 0 � 2 × 10 2 λ, this broadening effect reduces, and the behavior of angular distributions is more similar to the behavior in wave zone.

Oblique
Incidence. e angular distributions become asymmetric when a bunch with oblique incidence is considered. In wave zone, the asymmetry is more severe as the transverse beam size is large. In the prewave zone, however, the location in center of detector does have value (as shown in Figures 8 and 9); the value becomes large when the transverse beam size becomes large. e above results can be understood from a simple physical diagram as shown in Figure 10.

Case Study.
Although people normally consider a beam with incident angle of 45 degrees, in practice, however, in this study we consider the angular distribution from different incident angles. When a beam with incident angle is     considered, the angular distributions become asymmetric, as reported in Table 1 below. For 1 THz case, we can see that the asymmetry of angular distribution becomes severe as the incident angle is increased (as shown in Figure 11). For 10 THz case, when the incident angle is increased, the asymmetry pattern changes. e cause of such asymmetric angular distributions can be understood from a simple diagram as depicted in Figure 12.
In the above figure, we assume that θ 2 and θ 1 are the position with peak value, where θ 2 > 0 and θ 2 < 0. We can easily see the phase difference of the radiations from p 1 and p 2 is different at θ 2 and θ 1 . is causes asymmetric angular distributions.

Conclusions and Future Work
In this study, we proposed a theory treating prewave zone OTR with the consideration of beam structure. In the proposed theory, the beam transverse size effect can be simplified to transverse form factor when wave zone approximation is considered. A parallelized numerical integration code based on quasi-Monte Carlo method was developed to perform the complicated integration. We demonstrated, theoretically, that the angular distributions in the prewave zone are more sensitive than those in the wave zone. e simulations also provide consistent findings and results. e provided formalism theory and numerical tool could be useful for the prediction of radiation properties for the OTR experiment in prewave zone. e disadvantage of this method is that a large amount of samplings may need to be employed to achieve good convergence. In the future, we will investigate more robust technique to improve the performance and convergence speed using various parallelism concepts like the one used in this paper, i.e., MPI. Similarly, statistical methods should be used to find accurate mapping and excluding similar samples that could potentially decrease the convergence speed.
Data Availability e data can be requested from the corresponding author.

Conflicts of Interest
e authors declare that they have no conflicts of interest.