#### Abstract

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 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. 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.

#### 1. Introduction

When an electron traverses across an interface of two different media, radiation is emitted. This phenomenon is called transition radiation (TR) and was theoretically predicted by Ginzburg and Frank in 1946. Their approach to tackle this problem is based on solving the electromagnetic wave equations directly with consideration of boundary conditions at the interface [1–3] and provided solutions in wave zone (far-field). 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–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 fields 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. Kirschhoff’s method can be used to solve this scattering problem [8, 9]. With this integral formalism, one may find that the condition for far-field approximation holds when the observation distance *r*_{0} >> *γ*^{2}*λ*. In this regard, the radiation source is the surface current on conductor which is induced by the electric field from incident electron; the size, or more specifically the cut-off size, is of the order *γ*^{λ} because the asymptotic behavior of particle field falls off with radial distance *d* perpendicular to the trajectory that is proportional to exp(*−*2*πd/βγλ*). Due to the finite source size, the prewave zone can be defined as *γ*^{2}*λ* > *r*_{0} >> *γλ.* The behaviors of the OTR in this region, such as spectral angular distribution, are quite different from those in wave zone [10, 11].

The 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–18]. For a bunch with Gaussian distribution , the OTR angular spectrum is proportional to , where *I*_{e} is the OTR angular spectrum of single particle with normal incidence. Moreover, and 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. The 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. The 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. The 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). The 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) The developed Fortran program is then used to tackle the complex numerical integration problem. (iv) We parallelize the program with MPI programming interface.

The 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.

#### 2. 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. The beam size effects in prewave zone and wave zone are also discussed [19].

##### 2.1. Kirchhoff’s Method for OTR

In this section, we would review the theoretical treatment proposed by Karlovets [20]. The coordinate system in original work is in left-hand coordinate system and we will follow this convention.

###### 2.1.1. 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 are smooth function, *V* is a volume in space, *S* is the boundary enclosing *V*, and 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 is the mirror image of with respect to plane at *z* = 0. If we put we can get the fields relation for each component of **E**^{R}. To apply the integral equation to scattering problem, we can make and this yields where **E** and **E**_{0} are total field and incident field. With the boundary conditions of perfect conducting plane, , we can further conclude

**(a)**

**(b)**

For *r*_{0}>> *r,* we have and therefore the radiation fields (scattered field) become as illustrated by

###### 2.1.2. 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. Therefore, one should put the center of detector (COD) with angle *α* relative to *z*-axis. The 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 byand 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:

##### 2.2. 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. The moving fields from the charge particle are scattered by the media. This concept was proposed by Fermi [21]. When tackling a scattering problem, one normally considers the incident fields in frequency domain. Therefore, 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 having the distribution , where *N* is the number of electrons in this bunch, and 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:

The fields in the frequency domain arewhere is the longitudinal form factor for this bunch. The 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 isand velocity fields from this electron bunch in the target’s frame are

##### 2.3. 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. This yields the radiation fields in the target’s frame

Transforming to observer’s frame by using equation (9), we get

The following formula can be calculated by equations (11) and (16):

Substituting , we get

With the radiation fields, one can calculate the radiation energy with (10).

##### 2.4. 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 and ; 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. The 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. This is a time-consuming task, and therefore we parallelize this implementation code with the MPI programming interface and model. This will reduce the time needed to reach particular decisions. The detailed description for the algorithm and source code are given later in the appendices.

#### 3. 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.

##### 3.1. Finite Difference Time Domain Method

The 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).

**(a)**

**(b)**

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. This is the so-called perfectly matched layer. In the FDTD code, the uniaxial perfectly matched layer (UPML) is implemented. The 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. The region 2 is full of anisotropic medium which is characterized by

The 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 regionswhere Γ and *T* are reflection and transmission coefficient. By using the tangential boundary conditionswe can solve the reflection and transmission coefficient and also get Snell’s law

If *k*_{2x} *=* *bk*_{1x}

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 aswhere is

For two-dimension case, we have three kinds of material properties for UPML (Figure 3): (i) Boundary normal to *x* direction (*R*1): set *σ*_{y} = *σ*_{z} = 0. (ii) Boundary normal to *y* direction (*R*1): set *σ*_{x} = *σ*_{z} = 0. (iii) Overlapping region (*R*3): set *σ*_{z} = 0.

##### 3.2. 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. The code we use in this study is written in Fortran with MPI parallelization and is developed by Dr. Toseo Moritaka. The 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. The 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. The absolute square of the fields from different angles in frequency domain is corresponding to the angular distribution of radiation.

#### 4. Results and Discussion

##### 4.1. Theoretical Calculations

In this section, we show the radiation angular distribution from electron bunch with normal incidence. We also compare the angular distributions in wave zone with the angular distributions in prewave zone.

###### 4.1.1. Wave Zone vs. Prewave Zone

We studied the angular distributions from various transverse beam sizes at different observation distances. The 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. This shows that the angular distributions in prewave zone are more sensitive to transverse beam size.

**(a)**

**(b)**

##### 4.2. Simulations

To verify the theory, few two-dimensional FDTD simulations were performed. This 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).

###### 4.2.1. Evolution of *B*_{z} Field in Full Simulation Domain

###### 4.2.2. Simulations vs. Theory

To verify the theory, we performed the simulation for bunch with normal incidence and compared it with the results from theory. The obtained results and findings are shown in Figure 7.

**(a)**

**(b)**

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 *r*_{0} *=* 2 × 10^{2}*λ*, this broadening effect reduces, and the behavior of angular distributions is more similar to the behavior in wave zone.

##### 4.3. Oblique Incidence

The 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.

The above results can be understood from a simple physical diagram as shown in Figure 10.

##### 4.4. 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. The 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}. This causes asymmetric angular distributions.

#### 5. 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. The simulations also provide consistent findings and results. The provided formalism theory and numerical tool could be useful for the prediction of radiation properties for the OTR experiment in prewave zone. The 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

The data can be requested from the corresponding author.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the National Natural Science Foundation of China, 41774140, project title: Research on Soil Radon Measurement Technology and Abnormal Information Extraction Method in Uranium Exploration; National Natural Science Foundation of China, 11675028, project title: Research on Digital Shaping Algorithm and Real-Time Processing Technology of Nuclear Pulse Signal; and National Key R&D Program, 2018YFC0603304, project title: Development of Neutron Gamma Logging Methods and Instruments for Rock-Forming Elements.