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

Authors: Haofeng Liu; Fang Fang
Journal: Scientific Programming
 (2022)

---

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


---

## Body

## 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 r0 >> γ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λ > r0 >> γλ. 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 exp−x2+y2/2σ⊥2exp−z2/2σ∥2, the OTR angular spectrum is proportional to ℱ⊥λ,θℱ∥λIeλ,θA, where Ie 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 Section2, 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 (Figure1), 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.(1)∫Vf∇2g−g∇2fdV=∮Sn^⋅f∇g−g∇fdS,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:(2)g=−expikr−r0r−r0+expikr−r0′r−r0′,where r0′ is the mirror image of r0 with respect to plane at z = 0. If we put f=EiR we can get the fields relation for each component of ER. To apply the integral equation to scattering problem, we can make ERr0,ω=Er0,ω−E0r0,ω and this yields (3)ERr0,ω=14π∫E0r,ωn^⋅∇gdS,where E and E0 are total field and incident field. With the boundary conditions of perfect conducting plane, Ex,yr,ω|S=0, we can further conclude (4)Ex,yRr0,ω=−14π∫Ex,y0r,ωn^⋅∇gdS.Figure 1
Projections of the coordinate system on y-z and x-z plans. (a) Projection ony-z plan. (b) Projection on x-z plan.
(a)(b)Forr0>> r, we have (5)n^⋅∇g|S≈2ikz0expikr0−rr0−r,and therefore the radiation fields (scattered field) become as illustrated by(6)Ex,yRr0,ω=−2ikz04π∫Ex,y0r,ωexpikr0−rr0−rdS.


#### 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 by(7)ℛy⋅ℛx=cosθx0sinθx010−sinθx0cosθx1000cosα+θy−sinα+θy0−sinα+θycosα+θy,and the radiation fields transformation between these two frames is illustrated as(8)EkR=cosθxsinα+θysinθxcosα+θysinθx0cosα+θysinα+θy−sinθx−sinα+θycosθxcosα+θycosθx⋅ER.Assuming that the longitudinal component in the observer’s frame is absent, we get(9)ExkR=ExR,EykR=EyRcosα+θy+ExRtanθxtanα+θy.With the fields in observer’s frame, the radiation energy spectrum can be calculated with the following equation:(10)WdωdΩ=cr02Exk2+Eyk2.


### 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 + ze direction, we have(11)Exee=eωπv2γxexe2+ye2K1ωxe2+ye2vγexpiωzev,Eyee=eωπv2γyexe2+ye2K1ωxe2+ye2vγexpiωzev,Ezee=−eωπv2γiγK0ωxe2+ye2vγexpiωzev.When an electron moves inze direction in space with constant velocity, then there exists a time t0, such that this electron would arrive the plane ze = 0 when t = t0. If this electron has a deviation in transverse coordinate, say for example (x0, y0), the velocity fields in the time domain can be written as given by(12)Exe,ye,zeexe−x0,ye−y0,ze,t−t0.For a bunch of electrons moving alongze direction with constant velocity v having the distribution NB⊥x′,y′B∥t′, 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:(13)Exe,ye,zeb=N∭B⊥x′,y′B∥t′Exe,ye,zeexe−x′,ye−y′,ze,t−t′dx′dy′dt′.The fields in the frequency domain are(14)Exe,ye,zeb=Nℱ∥∫∫B⊥x′,y′Exe,ye,zeexe−x′,ye−y′,ze,ωdx′dy′,where ℱ∥=∫B∥t′expiωt′dt′ 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 is(15)xe=x,ye=ycosα−zsinα,ze=ysinα+zcosα,Exb=Exeb,Eyb=Eyebcosα+Ezebsinα,Ezb=−Ezebsinα+Ezebcosα,and velocity fields from this electron bunch in the target’s frame are(16)Exbxe,ye,ze,ω=Exebx,ycosα−zsinα,zcosα+ysinα,ω,Eybxe,ye,ze,ω=Eyebx,ycosα−zsinα,zcosα+ysinα,ωcosαn+Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinα,Ezbxe,ye,ze,ω=−Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinα+Ezebx,ycosα−zsinα,zcosα+ysinα,ωcosα.


### 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 (17)ExRr0,ω=−2ikz04π∫Exebx,ycosα−zsinα,zcosα+ysinα,ωexpikr0−rr0−rdS,EyRr0,ω=−2ikz04π∫Eyebx,ycosα−zsinα,zcosα+ysinα,ωcosα+Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinαexpikr0−rr0−rd.Transforming to observer’s frame by using equation (9), we get(18)ExkR=ExR=−2ikz04π∫Exebx,y,z,ωexpikr0−rr0−rdSEykR=EyRcosα+θy+ExRtanθxtanα+θy,=−2ikz04π∫1cosα+θyEyebx,ycosα−zsinα,zcosα+ysinα,ωcosα+Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinα+Exebx,y,z,ωtanθxtanα+θyexpikr0−rr0−rdS.The following formula can be calculated by equations (11) and (16):Substitutingx¯=x−x′,y¯=y−y′/cosα, we get(19)ExkR=Nℱ∥2ieβγλ2cosα+θyr0∭∫b⊥x′,y′x¯x¯2+y¯2cos2αK12πβγλx¯2+y¯2cos2αexpiϕdx¯dy¯dx′dy′,EykR=Nℱ∥2ieβγλ2cosθxr0∭∫b⊥x′,y′K12πβγλx¯2+y¯2cos2αy¯cos2α+x¯tanθxsinα+θyx¯2+y¯2cos2α−iγK02πβγλx¯2+y¯2cos2αsinαexpiϕdx¯dy¯dx′dy′,ϕ=kx¯+x′−x02+y¯+y′cosα−y02+z02+ky¯+y′/cosαsinαβ.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 ρ¯2=x¯2+y¯2cos2α and ρ′=x′2+y′2; this yields (20)ExkR=ℱ∥2ieβγλ212σ⊥21r0cosα+θycosα∭∫ρ¯ρ′exp−ρ′22σ⊥2K12πβγλρ¯expikΦdρ¯dρ′dϕ¯dϕ,′EykR=ℱ∥2ieβγλ212σ⊥21r0cosθxcosα∭∫ρ¯ρ′exp−ρ′22σ⊥2K12πβγλρ¯cosαsinϕ¯−iγK02πβγλρ¯sinαexpikΦdρ¯dρ′dϕ¯dϕ,′Φ=r0−1cosαρ¯sinϕ¯+ρ′sinϕ′sinα+θy−sinαβ+12r0ρ¯cosϕ¯+ρ′cosϕ′2+1cos2αρ¯sinϕ¯+ρ′sinϕ′2cos2α+θy.ℱ∥=1/2σ∥22πσ∥2exp−σ∥2ω2/2.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.


## 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 (Figure1), 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.(1)∫Vf∇2g−g∇2fdV=∮Sn^⋅f∇g−g∇fdS,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:(2)g=−expikr−r0r−r0+expikr−r0′r−r0′,where r0′ is the mirror image of r0 with respect to plane at z = 0. If we put f=EiR we can get the fields relation for each component of ER. To apply the integral equation to scattering problem, we can make ERr0,ω=Er0,ω−E0r0,ω and this yields (3)ERr0,ω=14π∫E0r,ωn^⋅∇gdS,where E and E0 are total field and incident field. With the boundary conditions of perfect conducting plane, Ex,yr,ω|S=0, we can further conclude (4)Ex,yRr0,ω=−14π∫Ex,y0r,ωn^⋅∇gdS.Figure 1
Projections of the coordinate system on y-z and x-z plans. (a) Projection ony-z plan. (b) Projection on x-z plan.
(a)(b)Forr0>> r, we have (5)n^⋅∇g|S≈2ikz0expikr0−rr0−r,and therefore the radiation fields (scattered field) become as illustrated by(6)Ex,yRr0,ω=−2ikz04π∫Ex,y0r,ωexpikr0−rr0−rdS.


### 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 by(7)ℛy⋅ℛx=cosθx0sinθx010−sinθx0cosθx1000cosα+θy−sinα+θy0−sinα+θycosα+θy,and the radiation fields transformation between these two frames is illustrated as(8)EkR=cosθxsinα+θysinθxcosα+θysinθx0cosα+θysinα+θy−sinθx−sinα+θycosθxcosα+θycosθx⋅ER.Assuming that the longitudinal component in the observer’s frame is absent, we get(9)ExkR=ExR,EykR=EyRcosα+θy+ExRtanθxtanα+θy.With the fields in observer’s frame, the radiation energy spectrum can be calculated with the following equation:(10)WdωdΩ=cr02Exk2+Eyk2.


## 2.1.1. Kirchhoff’s Method for Plane Conductor

Consider an OTR system (Figure1), 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.(1)∫Vf∇2g−g∇2fdV=∮Sn^⋅f∇g−g∇fdS,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:(2)g=−expikr−r0r−r0+expikr−r0′r−r0′,where r0′ is the mirror image of r0 with respect to plane at z = 0. If we put f=EiR we can get the fields relation for each component of ER. To apply the integral equation to scattering problem, we can make ERr0,ω=Er0,ω−E0r0,ω and this yields (3)ERr0,ω=14π∫E0r,ωn^⋅∇gdS,where E and E0 are total field and incident field. With the boundary conditions of perfect conducting plane, Ex,yr,ω|S=0, we can further conclude (4)Ex,yRr0,ω=−14π∫Ex,y0r,ωn^⋅∇gdS.Figure 1
Projections of the coordinate system on y-z and x-z plans. (a) Projection ony-z plan. (b) Projection on x-z plan.
(a)(b)Forr0>> r, we have (5)n^⋅∇g|S≈2ikz0expikr0−rr0−r,and therefore the radiation fields (scattered field) become as illustrated by(6)Ex,yRr0,ω=−2ikz04π∫Ex,y0r,ωexpikr0−rr0−rdS.


## 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 by(7)ℛy⋅ℛx=cosθx0sinθx010−sinθx0cosθx1000cosα+θy−sinα+θy0−sinα+θycosα+θy,and the radiation fields transformation between these two frames is illustrated as(8)EkR=cosθxsinα+θysinθxcosα+θysinθx0cosα+θysinα+θy−sinθx−sinα+θycosθxcosα+θycosθx⋅ER.Assuming that the longitudinal component in the observer’s frame is absent, we get(9)ExkR=ExR,EykR=EyRcosα+θy+ExRtanθxtanα+θy.With the fields in observer’s frame, the radiation energy spectrum can be calculated with the following equation:(10)WdωdΩ=cr02Exk2+Eyk2.


## 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 + ze direction, we have(11)Exee=eωπv2γxexe2+ye2K1ωxe2+ye2vγexpiωzev,Eyee=eωπv2γyexe2+ye2K1ωxe2+ye2vγexpiωzev,Ezee=−eωπv2γiγK0ωxe2+ye2vγexpiωzev.When an electron moves inze direction in space with constant velocity, then there exists a time t0, such that this electron would arrive the plane ze = 0 when t = t0. If this electron has a deviation in transverse coordinate, say for example (x0, y0), the velocity fields in the time domain can be written as given by(12)Exe,ye,zeexe−x0,ye−y0,ze,t−t0.For a bunch of electrons moving alongze direction with constant velocity v having the distribution NB⊥x′,y′B∥t′, 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:(13)Exe,ye,zeb=N∭B⊥x′,y′B∥t′Exe,ye,zeexe−x′,ye−y′,ze,t−t′dx′dy′dt′.The fields in the frequency domain are(14)Exe,ye,zeb=Nℱ∥∫∫B⊥x′,y′Exe,ye,zeexe−x′,ye−y′,ze,ωdx′dy′,where ℱ∥=∫B∥t′expiωt′dt′ 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 is(15)xe=x,ye=ycosα−zsinα,ze=ysinα+zcosα,Exb=Exeb,Eyb=Eyebcosα+Ezebsinα,Ezb=−Ezebsinα+Ezebcosα,and velocity fields from this electron bunch in the target’s frame are(16)Exbxe,ye,ze,ω=Exebx,ycosα−zsinα,zcosα+ysinα,ω,Eybxe,ye,ze,ω=Eyebx,ycosα−zsinα,zcosα+ysinα,ωcosαn+Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinα,Ezbxe,ye,ze,ω=−Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinα+Ezebx,ycosα−zsinα,zcosα+ysinα,ωcosα.


## 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 (17)ExRr0,ω=−2ikz04π∫Exebx,ycosα−zsinα,zcosα+ysinα,ωexpikr0−rr0−rdS,EyRr0,ω=−2ikz04π∫Eyebx,ycosα−zsinα,zcosα+ysinα,ωcosα+Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinαexpikr0−rr0−rd.Transforming to observer’s frame by using equation (9), we get(18)ExkR=ExR=−2ikz04π∫Exebx,y,z,ωexpikr0−rr0−rdSEykR=EyRcosα+θy+ExRtanθxtanα+θy,=−2ikz04π∫1cosα+θyEyebx,ycosα−zsinα,zcosα+ysinα,ωcosα+Ezebx,ycosα−zsinα,zcosα+ysinα,ωsinα+Exebx,y,z,ωtanθxtanα+θyexpikr0−rr0−rdS.The following formula can be calculated by equations (11) and (16):Substitutingx¯=x−x′,y¯=y−y′/cosα, we get(19)ExkR=Nℱ∥2ieβγλ2cosα+θyr0∭∫b⊥x′,y′x¯x¯2+y¯2cos2αK12πβγλx¯2+y¯2cos2αexpiϕdx¯dy¯dx′dy′,EykR=Nℱ∥2ieβγλ2cosθxr0∭∫b⊥x′,y′K12πβγλx¯2+y¯2cos2αy¯cos2α+x¯tanθxsinα+θyx¯2+y¯2cos2α−iγK02πβγλx¯2+y¯2cos2αsinαexpiϕdx¯dy¯dx′dy′,ϕ=kx¯+x′−x02+y¯+y′cosα−y02+z02+ky¯+y′/cosαsinαβ.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 ρ¯2=x¯2+y¯2cos2α and ρ′=x′2+y′2; this yields (20)ExkR=ℱ∥2ieβγλ212σ⊥21r0cosα+θycosα∭∫ρ¯ρ′exp−ρ′22σ⊥2K12πβγλρ¯expikΦdρ¯dρ′dϕ¯dϕ,′EykR=ℱ∥2ieβγλ212σ⊥21r0cosθxcosα∭∫ρ¯ρ′exp−ρ′22σ⊥2K12πβγλρ¯cosαsinϕ¯−iγK02πβγλρ¯sinαexpikΦdρ¯dρ′dϕ¯dϕ,′Φ=r0−1cosαρ¯sinϕ¯+ρ′sinϕ′sinα+θy−sinαβ+12r0ρ¯cosϕ¯+ρ′cosϕ′2+1cos2αρ¯sinϕ¯+ρ′sinϕ′2cos2α+θy.ℱ∥=1/2σ∥22πσ∥2exp−σ∥2ω2/2.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. (21)∇×E+∂B∂t=0,∂D∂t−∇×H=J,B=μH,D=ϵE.For two-dimensional system with the coordinate (x, y), the fields behavior can be split into transverse electric (TEz) mode (22)Hx=Hy=0,Ez=0,−μ∂Hz∂t=∂Ex∂y−∂Ey∂xϵ∂Ex∂t∂Hz∂y−Jx,ϵ∂Ey∂t=−∂Hz∂x−Jyn,and transverse magnetic (TMz) mode(23)Ex=Ey=0,Hz=0,ϵ∂Ez∂t=∂Hy∂x−∂Hx∂y−Jz,μ∂Hx∂t=−∂Ez∂y,μ∂Hy∂t=∂Ez∂x.One may then use central difference scheme to discretize Maxwell’s equation and get the following discretization for TEz mode and for TMz mode (Figure 2).Figure 2
Two-dimension Yee grid for both TE and TM modes. (a) TE mode for 2D Yee grid. (b) TM mode for 2D Yee grid.
(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 (24)Hinc=H0exp−jk1xx−jk1yy,which propagates toward region 2. The region 2 is full of anisotropic medium which is characterized by(25)ϵ2↔=ϵ2a000b000b,μ2↔=μ2c000d000d.The Maxwell equations in this medium are(26)∇×E=−jωμ2↔H,∇×H=jωϵ2↔E.We consider here the case forTEz. To calculate the reflection coefficient, we write down the plane wave in both regions(27)H1=z^H01+Γexp2jk1xxexp−jk1xx−jk1yy,E1=H0−x^k1yωϵ11+Γexp2jk1xx+y^k1xωϵ11−Γexp2jk1xxexp−jk1xx−jk1yy,H2=z^H0Texp−jk2xx−jk2yy,E2=H0−x^k2yωϵ2a+y^k2xωϵ2bTexp−jk1xx−jk1yy,where Γ and T are reflection and transmission coefficient. By using the tangential boundary conditions(28)H1∥=H2∥,E1∥=E2∥,we can solve the reflection and transmission coefficient (29)Γ=k1x−k2xb−1k1x+k2xb−1,and also get Snell’s law (30)k2x=k22bd−k1y2a−1b.Ifϵ1=ϵ2,μ1=μ2,d=banda−1=b,thenk1=k2andk2x = bk1xSubstituting 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 tox 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(31)∇×E=−jωμs↔H,∇×H=jωϵstE,where s↔ is(32)st=sx−1sysz000sxsy−1sz000sxsysz−1.For two-dimension case, we have three kinds of material properties for UPML (Figure3):(i) Boundary normal tox direction (R1): set σy = σz = 0.(ii) Boundary normal toy direction (R1): set σx = σz = 0.(iii) Overlapping region (R3): set σz = 0.Figure 3
UPML in two-dimension simulation.


### 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 Figure4.Figure 4
Simulation setup for OTR from normal incidence of electron bunch.Once the current pulse moves across the PEC boundary, OTR will be generated and then propagates toward the right hand side. The time signals forHz 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.


## 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. (21)∇×E+∂B∂t=0,∂D∂t−∇×H=J,B=μH,D=ϵE.For two-dimensional system with the coordinate (x, y), the fields behavior can be split into transverse electric (TEz) mode (22)Hx=Hy=0,Ez=0,−μ∂Hz∂t=∂Ex∂y−∂Ey∂xϵ∂Ex∂t∂Hz∂y−Jx,ϵ∂Ey∂t=−∂Hz∂x−Jyn,and transverse magnetic (TMz) mode(23)Ex=Ey=0,Hz=0,ϵ∂Ez∂t=∂Hy∂x−∂Hx∂y−Jz,μ∂Hx∂t=−∂Ez∂y,μ∂Hy∂t=∂Ez∂x.One may then use central difference scheme to discretize Maxwell’s equation and get the following discretization for TEz mode and for TMz mode (Figure 2).Figure 2
Two-dimension Yee grid for both TE and TM modes. (a) TE mode for 2D Yee grid. (b) TM mode for 2D Yee grid.
(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 (24)Hinc=H0exp−jk1xx−jk1yy,which propagates toward region 2. The region 2 is full of anisotropic medium which is characterized by(25)ϵ2↔=ϵ2a000b000b,μ2↔=μ2c000d000d.The Maxwell equations in this medium are(26)∇×E=−jωμ2↔H,∇×H=jωϵ2↔E.We consider here the case forTEz. To calculate the reflection coefficient, we write down the plane wave in both regions(27)H1=z^H01+Γexp2jk1xxexp−jk1xx−jk1yy,E1=H0−x^k1yωϵ11+Γexp2jk1xx+y^k1xωϵ11−Γexp2jk1xxexp−jk1xx−jk1yy,H2=z^H0Texp−jk2xx−jk2yy,E2=H0−x^k2yωϵ2a+y^k2xωϵ2bTexp−jk1xx−jk1yy,where Γ and T are reflection and transmission coefficient. By using the tangential boundary conditions(28)H1∥=H2∥,E1∥=E2∥,we can solve the reflection and transmission coefficient (29)Γ=k1x−k2xb−1k1x+k2xb−1,and also get Snell’s law (30)k2x=k22bd−k1y2a−1b.Ifϵ1=ϵ2,μ1=μ2,d=banda−1=b,thenk1=k2andk2x = bk1xSubstituting 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 tox 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(31)∇×E=−jωμs↔H,∇×H=jωϵstE,where s↔ is(32)st=sx−1sysz000sxsy−1sz000sxsysz−1.For two-dimension case, we have three kinds of material properties for UPML (Figure3):(i) Boundary normal tox direction (R1): set σy = σz = 0.(ii) Boundary normal toy direction (R1): set σx = σz = 0.(iii) Overlapping region (R3): set σz = 0.Figure 3
UPML in two-dimension simulation.


## 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 Figure4.Figure 4
Simulation setup for OTR from normal incidence of electron bunch.Once the current pulse moves across the PEC boundary, OTR will be generated and then propagates toward the right hand side. The time signals forHz 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 Figure5. We can firstly observe that the angular distributions change significantly with distance and become constant when the distance is large enough, such as r0 = 5 × 103λ and r0 = 1 × 104λ. 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 r0 = 1 × 102λ, 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.Figure 5
Angular distributions for various beam transverse sizes at different observation distances. (b) Peak positions for profiles from (a).
(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 (Figure6).Figure 6
Heat maps ofBz field in full simulation domain.


#### 4.2.1. Evolution ofBz 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 Figure7.Figure 7
Comparison of angular spectrums between simulation and theory. (b) Peak positions for profiles from (a).
(a)(b)We can see the results from simulation qualitatively match the prediction form theory: in the prewave zone, sayr0 = 1 × 102λ, the angular distribution from larger beam transverse size can be broader than that from smaller beam transverse size. Also, as the distance increases, like r0 = 2 × 102λ, 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 Figures8 and 9); the value becomes large when the transverse beam size becomes large.Figure 8
Angular distributions withα = π/4 in wave zone.Figure 9
Angular distributions withα = π/4 in prewave zone.The above results can be understood from a simple physical diagram as shown in Figure10.Figure 10
A simple picture for OTR from electron bunch.


### 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 Table1 below.Table 1
Parameters from HBI group, NSRRC.
ProjectNumerical valueEnergy30 MevTransverse beam size (rms)2 mmObservation distance1 mFor 1 THz case, we can see that the asymmetry of angular distribution becomes severe as the incident angle is increased (as shown in Figure11). 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.Figure 11
Angular distributions with different observation angles.Figure 12
Simple picture of asymmetric angular distributions.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 p1 and p2 is different at θ2 and θ1. This causes asymmetric angular distributions.


## 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 Figure5. We can firstly observe that the angular distributions change significantly with distance and become constant when the distance is large enough, such as r0 = 5 × 103λ and r0 = 1 × 104λ. 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 r0 = 1 × 102λ, 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.Figure 5
Angular distributions for various beam transverse sizes at different observation distances. (b) Peak positions for profiles from (a).
(a)(b)


## 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 Figure5. We can firstly observe that the angular distributions change significantly with distance and become constant when the distance is large enough, such as r0 = 5 × 103λ and r0 = 1 × 104λ. 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 r0 = 1 × 102λ, 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.Figure 5
Angular distributions for various beam transverse sizes at different observation distances. (b) Peak positions for profiles from (a).
(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 (Figure6).Figure 6
Heat maps ofBz field in full simulation domain.


### 4.2.1. Evolution ofBz 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 Figure7.Figure 7
Comparison of angular spectrums between simulation and theory. (b) Peak positions for profiles from (a).
(a)(b)We can see the results from simulation qualitatively match the prediction form theory: in the prewave zone, sayr0 = 1 × 102λ, the angular distribution from larger beam transverse size can be broader than that from smaller beam transverse size. Also, as the distance increases, like r0 = 2 × 102λ, this broadening effect reduces, and the behavior of angular distributions is more similar to the behavior in wave zone.


## 4.2.1. Evolution ofBz 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 Figure7.Figure 7
Comparison of angular spectrums between simulation and theory. (b) Peak positions for profiles from (a).
(a)(b)We can see the results from simulation qualitatively match the prediction form theory: in the prewave zone, sayr0 = 1 × 102λ, the angular distribution from larger beam transverse size can be broader than that from smaller beam transverse size. Also, as the distance increases, like r0 = 2 × 102λ, 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 Figures8 and 9); the value becomes large when the transverse beam size becomes large.Figure 8
Angular distributions withα = π/4 in wave zone.Figure 9
Angular distributions withα = π/4 in prewave zone.The above results can be understood from a simple physical diagram as shown in Figure10.Figure 10
A simple picture for OTR from electron bunch.


## 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 Table1 below.Table 1
Parameters from HBI group, NSRRC.
ProjectNumerical valueEnergy30 MevTransverse beam size (rms)2 mmObservation distance1 mFor 1 THz case, we can see that the asymmetry of angular distribution becomes severe as the incident angle is increased (as shown in Figure11). 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.Figure 11
Angular distributions with different observation angles.Figure 12
Simple picture of asymmetric angular distributions.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 p1 and p2 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.



---

*Source: 6480380-2022-05-30.xml*