[label1]Shandong Institute of Advanced Technology, 1501 Panlong Road, Jinan
250100, China [label2]Institute of High Energy Physics, Chinese Academy of
Sciences, 19B Yuquan Road, Shijingshan District, Beijing 100049, China
[label3]School of Physics Sciences, Institute of Physics, University of
Chinese Academy of Sciences, Chinese Academy of Sciences, 3rd South
Zhongguancun Street, Haidian District, Beijing 100190, China [label4]Jilin
University, 2699 Qianjin Street, Changchun 130216, China

[label5]Beijing National Laboratory for Condensed Matter Physics, Institute of
Physics, University of Chinese Academy of Sciences, Chinese Academy of
Sciences, 3rd South Zhongguancun Street, Haidian District, Beijing 100190,
China

[label6]Liaoning University, 66 Chongshanzhong Road, Huanggu District,
Shenyang 110036, China [label7]Shandong University, 27 Shanda Nanlu, Jinan
250100, China

# Electric Field Measurement by Edge Transient Current Technique on Silicon
Low Gain Avalanche Detector

Chenxi Fu  Haobo Wang  Ryuta Kiuchi  Jianing Lin  Xin Shi  Tao Yang  Xiaoshen
Kang  Weimin Song  Congcong Wang  Weiwei Xu  Zijun Xu  Suyu Xiao
<EMAIL_ADDRESS>
###### Abstract

A novel methodology, named the diffusion profile method, is proposed in this
research to measure the electric field of a low gain avalanche detector
(LGAD). The proposed methodology utilizes the maximum of the time derivative
of the edge transient current technique (edge-TCT) test waveform to quantify
the dispersion of the light-induced carriers. This method introduces the
estimation of the elongation of the carrier cluster caused by diffusion and
the divergence of the electric field force during its drift along the
detector. The effectiveness of the diffusion profile method is demonstrated
through the analysis of both simulated and measured edge-TCT waveforms.
Experimental data was collected from a laser scan performed on an LGAD
detector along its thickness direction. A simulation procedure has been
developed in RASER (RAdiation SEmiconductoR) to generate signals from LGAD. An
assumption of immediate one-step carrier multiplication is introduced to
simplify the avalanche process. Simulation results were compared with
transient current data at the waveform level and showed a favorable match.
Both simulation and experimental results have shown that the diffusion profile
method could be applied to certain edge-TCT facilities as an alternative of
electric field measurement.

###### keywords:

low gain avalanche detector , edge transient current technique , diffusion
profile

##  1 Introduction

High-energy experiments, such as the High Luminosity Large Hadron Collider
(HL-LHC) [12, 13, 14], undergo periodic updates. In order to meet the
requirements of high luminosity facilities [15, 16], there has been a
significant increase in the need to enhance and improve particle trackers,
particularly in terms of time performance. Low Gain Avalanche Detectors (LGAD)
[17] have gained significant attention as potential candidates for various
applications due to their exceptional timing capabilities. It has been
reported that Silicon LGAD are able to achieve timing resolutions as low as 20
ps [18, 19]. In order to ensure the optimal performance of LGAD, it is crucial
to carry out comprehensive testing and measurements before their
implementation [19, 20].

The Transient Current Technique (TCT) is a method used to generate a signal on
particle detectors for the purpose of investigating their properties. This
technique has been employed for various purposes, such as determining the
effective doping level of detectors [29] and evaluating the gain rate of LGAD
[30]. Additionally, the determination of the built-in electric field in p-i-n
planar and strip detectors can be achieved by modifying the direction of the
injected laser using the Edge Transient Current Technique (Edge-TCT) [10], by
analyzing the carrier velocity profile [11]. Nevertheless, the direct
application of this method to LGAD is not feasible due to specific
limitations. Instead, this study introduces a novel approach called the
diffusion profile method. TCT generates non-equilibrium carrier clusters
within the detector using a laser. Before the carrier avalanche, the
distribution of clusters undergoes modifications during the drift process.
This phenomenon is the result of diffusion and the divergence of the electric
field force. A flattened carrier distribution results in a signal with a
weaker peak, indicating a lower intensity of the electric field or a higher
degree of field variation.

The RAdiation SEmiconductoR (RASER) is a Python package designed to simulate
physics processes in semiconductor detectors [raser]. The software provides
extensibility and customization, distinguishing it from other detector
simulators such as Allpix Squared [26], Weightfield2 [27], and TRACS [23],
which are coded in C++. The package has shown strong agreement with
experimental data concerning the temporal performance of non-irradiated p-i-n
diodes [6]. Furthermore, previous studies utilizing this package have shown
successful predictions for 3D detectors [7]. The aim of employing RASER
simulation in this research is to replicate the essential steps of the edge-
TCT process in p-i-n and LGAD in order to verify the effectiveness of the
diffusion profile method and establish a foundation for future research. This
study represents the initial attempt to model an LGAD as an active device
capable of responding to laser beam input. A simplified model for carrier
multiplication in an LGAD has been successfully developed in this study. Our
approach is based on the underlying assumption that a single carrier entering
the gain layer generates $M$ carrier pairs at a specific space-time
coordinate. Here, the variable $M$ represents the rate of gain.

##  2 Facilities

###  2.1 Detector under test

The LGAD being tested, along with the p-i-n detector used for comparison, was
manufactured by Hamamatsu Photonics K.K. (HPK) [2], specified as HPK Type
3.1-50. The p-type bulk of the detectors has an active thickness of $50\
\rm\mu m$ and a size of $1300\times 1300\ \rm\mu m^{2}$. Guard rings are
distributed outside the active area. The I-V and C-V characteristics were
measured under reverse bias, as shown in Figure 1(b). The experimental setup
consisted of a Keithley 2410 power supply, a Keithley 6487 current meter, and
an Agilent E4980A LCR meter operating at a frequency of $10\rm kHz$.

(a) I-V relation of p-i-n and LGAD. (b) C-V relation of p-i-n and LGAD.

Figure 1: I-V and C-V relations of p-i-n and LGAD under test.

The abrupt junction model has been used to estimate the distribution of
doping. Both the p-i-n and LGAD structures exhibit a bulk doping concentration
of $1\times 10^{12}\ \rm cm^{-3}$. Additionally, the LGAD structure
incorporates a deep-doped region with a concentration of $1.913\times 10^{18}\
\rm cm^{-3}$, which is estimated from C-V data and the gain rate. The deep-
doped area is specifically modeled within the LGAD structure, covering a range
of $1\sim 2\ \rm\mu m$.

###  2.2 Experiment Setup

Figure 2: Schematic view of the edge-TCT setup at CERN RD50 collaboration [3].
Experiment data recorded in 2019.

The experimental setup for the edge-TCT scan is illustrated in Figure 2. An
infrared laser with a wavelength of 1064 nm is utilized, which has an
attenuation length of approximately 1 mm in silicon. The given length is
sufficient for the laser beam to effectively penetrate the detector. The beam
waist of the pulse has a full width at half maximum (FWHM) measurement of 8
$\rm\mu$m, and the pulse duration is 350 ps. Denoting the direction of the
laser beam as the $x$-axis, the direction along the edge as the $y$-axis, and
the direction of the electric field inside the detector as the $z$-axis serves
as a convenient coordinate system for analysis. The value of z is initialized
to zero at the top of the detector. The laser beam should be positioned
perpendicular to the edge surface of the detector. Three portable tables have
been arranged to facilitate the adjustment of the laser’s focal point,
enabling a comprehensive z-scan across the entire thickness of the detector (0
to 50 $\rm\mu$m) with a step size of 1 $\rm\mu$m. Before capturing the laser-
induced waveform, the laser beam is precisely directed towards the edge of the
active area of the detector. This is achieved by adjusting the focus of the
laser in a direction perpendicular to the laser beam at different positions
along the z-axis. This ensures that the generated electron-hole pairs are
within a restricted range near the designated value of $z$. The resulting
variation in waveform amplitude is then observed to determine the minimum
width of the laser beam. A Keithley 2410 source is used as a high-voltage
power supply. To enhance the signal, a CIVIDEC [4] wideband amplifier with a
bandwidth of 2 GHz, a gain of 40 dB, and an integrated bias-Tee is utilized.
After shaping, the signal was converted into a digital format and recorded
using a Keysight DSO9254A oscilloscope [5], which has a bandwidth of 20 GHz
and a sample rate of 20 GSa/s.

In the following TCT tests, the devices were measured at room temperature.
Except for gain measurement, the devices were biased at 200 V from the
backside.

##  3 Modeling

The modeling of the edge-TCT process can be categorized into four main
components: field solution, generation of laser-induced non-equilibrium
carriers, collection of carrier drift signals, and filtering of the final
waveform. For LGAD, an additional step of carrier multiplication is
implemented.

###  3.1 Field profile

The electric potential $U$ of a detector can be calculated by solving the
Poisson equation

$\nabla^{2}U(\vec{x})=-eN_{eff}(\vec{x})/\varepsilon,$ (1)

in which $\vec{x}=(x,y,z)$ represents the vector coordinates, $e$ denotes the
unit charge, $\varepsilon$ denotes the dielectric constant of the detector
material, and $N_{eff}$ denotes the effective doping concentration. The
depletion ansatz is used to simplify the calculation. In addition, it is
assumed that the depletion area of the heavily-doped Ohmic contact layers is
thin enough to be considered negligible.

The calculation of the weighting potential for signal generation is performed
in a similar manner, using the following equation

$\nabla^{2}U_{w}(\vec{x})=0.$ (2)

At the readout electrode, $U_{w}$ is set to 1, while on all other electrodes,
$U_{w}$ is set to 0.

Within regions of high electric fields, carriers undergo acceleration and gain
enough energy to ionize electrons in the valence band. The multiplication
behavior of carriers is determined by the ionization rate, which is defined as
the number of electron-hole pairs generated by a drifting carrier as it
travels a unit distance. The ionization rate for electrons ($\alpha_{n}$) and
holes ($\alpha_{p}$), which is determined by the temperature $T$ and intensity
$E$, can be expressed as [8]

$\alpha_{p,n}(T,E)=\gamma(T)a_{p,n}\exp(-\frac{\gamma(T)b_{p,n}}{E}),$ (3)

where $a_{p,n}$ and $b_{p,n}$ are fixed parameters, and
$\gamma(T)=\tanh(\hbar\omega/2kT_{0})/\tanh\left[\hbar\omega/2kT(T/T_{0})\right]$
is a correction term for phonon-electron interaction, with $T_{0}=300\rm\ K$
and $\hbar\omega$ denoting the average phonon energy within the material [24].

###  3.2 Laser and carrier generation

The rate at which light-induced carriers are generated is proportional to the
intensity of the laser beam. For a typical Gaussian laser beam in the $x$
direction, the light intensity distribution is [25]

$I(r,x,t)=\frac{E_{p}}{\Delta t}\frac{4\sqrt{\ln
2}}{\pi^{3/2}w^{2}(x)}\exp{\frac{-2r^{2}}{w^{2}(x)}}\exp{\frac{-4t^{2}\ln
2}{(\Delta t)^{2}}}$ (4)

where $r^{2}=y^{2}+z^{2}$, $E_{p}$ and $\Delta t$ are the single pulse energy
and the pulse duration of the laser, and $w_{0}=\sqrt{\lambda/\pi
NA},w(x)=w_{0}\sqrt{1+(\lambda x/\pi w_{0}^{2}n)^{2}}$ are the beam widths at
$0$ and $x$ (with the origin at the beam focus), determined by the numerical
aperture ($NA$) of the focusing lens and the laser wavelength ($\lambda$). The
spatial FWHM of the laser beam is related to $w(z)$ by
$\mathrm{FWHM}_{spatial}=w(z)\sqrt{2\ln{2}}$.

Considering refraction and absorption, the intensity function will change its
form into

$I^{\prime}(r,x,t)=I_{0}(r,x,t)\exp\left[-\alpha(\lambda)x\right].$ (5)

The photons emitted by the laser beam are converted into excitation energy for
the carriers. During a small time interval $dt$, the density of generated
carriers can be calculated by

$\mathrm{d}N_{0}^{e-h\ pairs}=\alpha I^{\prime}\mathrm{d}t/h\nu,$ (6)

where $h\nu$ represents the energy of a photon.

###  3.3 Carrier drift and multiplication

Once the carriers are excited, they will drift through the bulk material,
thereby generating an electric current. According to the mobility formula [22]

$v_{p,n}=\mu_{p,n}(T,N_{eff},E)E,$ (7)

the velocity of a carrier at a given position can be mathematically
represented as a function of the effective doping concentration $N_{eff}$ and
the electric field intensity $E$. In this study, the Reggiani model [9] is
employed. Aside from oriented drifting, a Langevin term is considered to
simulate the random motion of the carriers by adding a Gaussian displacement
$\mathrm{d}\vec{x}_{Langevin}$, the probability distribution function of which
is

$f(\mathrm{d}\vec{x}_{Langevin})=\frac{1}{(2\pi
D\mathrm{d}t)^{3/2}}\exp{\left(-\frac{|\mathrm{d}\vec{x}_{Langevin}|^{2}}{2D\mathrm{d}t}\right)}$
(8)

where

$D=k_{B}T\mu/q$ (9)

is the Einstein diffusion coefficient. Here, $k_{B}$, $T$, and $q$ refer to
the Boltzmann constant, temperature, and elementary charge, respectively.

According to Shockley-Ramo’s theorem, a signal $i_{q}(t)$ will be generated at
the reading electrode during the drift of a charge $q$ by

$i_{q}(t)=q\vec{v}_{q}(t)\cdot\nabla U_{w}[\vec{x}_{q}(t)].$ (10)

The complete signal waveform can be acquired by combining the signals from all
generated carriers. Additionally, in LGAD, a carrier in the gain layer induces
carrier multiplication, which leads to signal amplification.

Instead of simulating the avalanche step by step, as illustrated in Figure 3,
it is assumed that the entry of a single carrier into the gain layer will
result in the generation of $M$ pairs of carriers without ionizing capability
at the bottom of the gain layer, where

$M=\frac{\exp\left[\int_{0}^{d_{gain}}\left(\alpha_{\mathrm{n}}-\alpha_{\mathrm{p}}\right)\mathrm{d}x\right]}{1-\int_{0}^{d_{gain}}\alpha_{\mathrm{p}}\exp\left[\int_{0}^{x}\left(\alpha_{\mathrm{n}}-\alpha_{\mathrm{p}}\right)\mathrm{d}x\right]\mathrm{d}x}$
(11)

is the gain factor of the device [22]. As the values of $\alpha_{n}$ and
$\alpha_{p}$ are sufficiently small outside the designated gain layer, it is
advantageous to only integrate within the gain layer, which has a width of
$d_{gain}$. The integral diverges when the device breaks down. For the LGAD
under test, the value of $M$ is determined to be 23.4.

Due to the double-exponential relationship between the field intensity $E$ and
the parameter $M$ (see Eq. 3 and 11), even a minor modification in the doping
profile can result in a substantial change in the gain rate. This statement
elucidates the challenges associated with the fabrication process of LGAD.
However, in simulation, it is possible to make a numerical adjustment to the
doping profile (by $10^{12}\rm\ cm^{-3}$) in order to align the gain rate
value with experimental results.

Figure 3: Sketch of simplified carrier multiplication process in RASER,
limited at the boundary between gain layer and the bulk. This picture is not
to scale.

###  3.4 Signal collection

Amplification is needed for the original induced currents to be detected,
stored, and analyzed. The electronic system can be modeled as shown in Figure
4, which incorporates the signal effects of various devices. These effects
include the capacitance $C_{D}$ and the AC resistance $R_{D}$ of the detector.
Additionally, the inductance $L_{T}$ and capacitance $C_{T}$ of the bias-tee
contribute to these effects. Furthermore, the gain ($A$) of the amplifier, the
input resistance ($R_{in}$), the output resistance ($R_{out}$), and the
capacitance of the cable ($C_{w}$) also contribute to these effects.

Figure 4: Equivalent circuit of Figure 2.

In the frequency domain, if $1/R_{D}$ and $1/L_{T}$ are small enough to be
neglected, the transfer function of the aforementioned devices can be
expressed as (letting $s=j\omega$)

$H(s)=\frac{R_{in}C_{T}}{C_{T}+C_{D}+R_{in}sC_{T}C_{D}}\frac{A}{sC_{w}R_{out}+1},$
(12)

and by using Laplace transformation, the amplification and shaping of the
original signal can be expressed through convolution with

$h(t)=A^{*}\frac{\exp(-t/\tau_{RC_{D}})-\exp(-t/\tau_{RC_{w}})}{\tau_{RC_{D}}-\tau_{RC_{w}}},$
(13)

where $\tau_{RC_{D}}=R_{in}C_{T}C_{D}/(C_{T}+C_{D})$ and
$\tau_{RC_{w}}=R_{out}C_{w}$ are the time constants of the circuit.
Additionally, $A^{*}=A\tau_{RC_{D}}/C_{D}$ is the total amplification rate,
measured in units of $\Omega$, that converts the original current signal into
voltage.

##  4 Waveform Simulation

This section aims to evaluate the results of the edge-TCT z-scan in the
direction of the device thickness. A global normalization factor is applied to
both p-i-n and LGAD. For LGAD, the laser intensity is 0.58 times stronger than
that of p-i-n in order to obtain a more distinct waveform.

In Figure 5, a good match of waveforms is exhibited when comparing the
simulation and experimental results of p-i-n and LGAD. But to ensure reliable
edge-TCT scan simulation, the waveform quantities, such as amplitude, rise
time, and charge collection efficiency, should have a good agreement with
experimental results at every point.

(a) Edge-TCT waveform of p-i-n. (b) Edge-TCT waveform of LGAD.

Figure 5: Edge-TCT waveform of both detectors, laser focus at $z$ direction of
$25\ \rm\mu m$, bias voltage 200 V.

(a) Amplitude-z relation.

(b) Rise time-z relation.

(c) Charge collection-z relation.

(d) Charge collection-Bias voltage relation.

Figure 6: Edge-TCT scan results of LGAD, top at $z=0$, bottom at $z=50\ \rm\mu
m$.

###  4.1 Amplitude

Figure 6(a) illustrates the relationship between the amplitude ($A$) of the
waveform and the position of the laser focus ($z$). The simulation results
demonstrate an overestimation in the range of $z\in(40,50)$ and an
underestimation in the range of $z\in(50,60)$. This phenomenon may potentially
be attributed to the diffusion of dopants inside the substrate, resulting in a
non-abrupt junction between the substrate and the bulk. Additionally, the lack
of consideration for light-induced carriers within the substrate, which
diffuse into the active layer, may have also played a role.

###  4.2 Rise time

The rise time of a signal can be determined by performing linear fitting on
the data points within the range of 20% to 80% of the signal’s amplitude. The
relationship between the rise time and the variable $z$ is illustrated in
Figure 6(b). An uniform result of 0.7 ns is observed and simulated for LGAD,
indicating stable signal generation and an excellent signal-to-noise ratio.
The presence of jaggedness within the figure can be attributed to the use of
discrete data points.

###  4.3 Charge collection

In this work, the charge collection of a waveform is defined as the integral
of the current over a specific time period. The result of the charge versus
focus $z$ coordinate is shown in Figure 6(c). As mentioned in Section 4.1,
apart from the junction and carrier diffusion at the bottom, the charge
collection successfully reflected the deposition of laser energy.

###  4.4 Gain measurement

To measure the gain rate of LGAD, we examined the relationship between charge
and reverse bias voltage $V$. Figure 6(d), which was sampled at $z=25\rm\ \mu
m$. No charge collection change is observed for p-i-n, thus the charge
collection of p-i-n is used as a $gain=1$ contract. Under the given voltage
range of $[60,200]\rm\ V$, the gain rate of LGAD increases from 4 to 26. At
200V, the integral in 11 is 23.4, while the simulation result from the
waveform is 24, with a relative error of less than 10%.

With comparisons of several quantities, an edge-TCT z-scan of silicon LGAD is
performed to ensure the accuracy of LGAD simulation. This work solidifies and
enhances the method of simplifying multiplication, making it more reliable. It
serves as a foundation for future research.

##  5 Electric Field Measurement

###  5.1 Velocity profile

The primary utility of edge-TCT lies in its ability to quantify the electric
field within the bulk of the detectors. The velocity profile method is widely
employed, particularly for p-i-n detectors [11]. For an electron or hole
moving within the interior of the detector’s active area, its velocity is
determined by the field intensity, as described by Equation 7. The signal,
which is primarily determined by the drifting velocity (excluding trap
capture), can be expressed using Equation 10. That is, it has the capability
to extract field intensity information from the edge-TCT waveform.

Considering a pair of carriers generated at a specific location $z=z_{0}$, the
sum of their velocities at the instant they are generated is

$\displaystyle v_{e}(t)+v_{h}(t)$
$\displaystyle=\mu_{e}[E(z_{e})]E(z_{e})+\mu_{h}[E(z_{h})]E(z_{h})$ (14)
$\displaystyle\approx\\{\mu_{e}[E(z_{0})]+\mu_{h}[E(z_{0})]\\}E(z_{0}).$

If a group of carriers conforms to the assumption that they are produced
within a limited region and time period around $t_{0}$, then the resulting
signal they generate should be

$i_{q}(z;t\to t_{0})=\sum_{q}q\vec{v}_{q}\nabla
U_{w}\approx\frac{ne}{d}(\mu_{e}+\mu_{h})E(z)$ (15)

as $\nabla U_{w}\approx(1/d)\ \hat{z}$, and $\vec{E}\approx E\hat{z}$ in
planar structure detectors. After amplification and RC shaping, it can be
observed that this proportional relationship may not hold strictly true.
However, the functional relationship $I(z;t\to t_{0})\approx I(E(z))$ should
remain valid.

To counteract the effect of carriers leaving the laser focus and carriers
emerging, the velocity profile is chosen as the signal value 0.1 ns after the
waveform emerges. The simulation results of the LGAD, as depicted in Figure 7,
demonstrate a significant peak at $z=5$ for the simulation and $z=10$ for the
experiment. Both results show no correlation with the known theoretical LGAD
field. The velocity profile method actually loses its effectiveness for LGAD
with our facilities.

(a) Velocity profile of LGAD.

(b) Diffusion profile of LGAD. The dashed box represents a preferred range
that excludes the effect of the non-Gaussian shape of the carrier cluster.

Figure 7: Two methods for LGAD electric field measurement. The pink region
represents the gain layer, while the cyan region represents the bulk.

This could actually be proven by the analyses provided below. According to
Equation 15, the extracted signal, denoted as $i_{q}$, should exclude the gain
signal. However, it is not possible to distinguish between the peak of laser
intensity and the peak of the gain signal (see Figure 5(b)). To provide a more
precise and comprehensive explanation, in systems where the time required for
carrier drifting is roughly equivalent to the duration of the laser pulse
(both of which are approximately 0.5 ns in our facilities), the ”velocity
profile” does not provide any information about the coordinates of the focal
point. Consequently, it is not suitable for quantifying velocity or field
intensity. Additionally, the back-end electronics may also impact the
resolution. For systems of this nature, a specific methodology is required to
analyze the internal built-in field, using additional variables. In addition
to considering the idealized signal input, an alternative approach involves
examining the collective behavior of non-equilibrium carriers induced by light
and analyzing their statistical behavior.

###  5.2 Diffusion Profile

In this section, the width of the propagating carrier cluster is being
researched. Ignoring carrier recombination, the distribution of a cluster of
light-induced non-equilibrium carriers drifting and diffusing under a uniform
electric field is [24]

$N=N_{0}\ast\frac{1}{\sqrt{2\pi\sigma_{D}^{2}}}\exp{[-\frac{(z-\mu
Et)^{2}}{2\sigma_{D}^{2}}]}$ (16)

where $N=N^{e-h\ pairs}(z,t)$, $N_{0}=N_{0}^{e-h\ pairs}(z,t)$ represents the
initial distribution of the carriers as defined in equation 5, and
$\sigma_{D}^{2}=2Dt$ if the diffusion coefficient $D$ is constant. The
convolution is performed in both the spatial and temporal dimensions.

For an initial spatial-temporal Gaussian distribution with spatial spread
$\sigma_{0}$ and temporal spread $\tau$, the convolution can be reduced to a
variance composition as

$\displaystyle N(z,t)$
$\displaystyle=\frac{1}{\sqrt{2\pi\tau^{2}}}\exp{(-\frac{t^{2}}{2\tau^{2}})}$
(17) $\displaystyle\ast\frac{\sum
N}{\sqrt{2\pi\sigma^{2}}}\exp{[-\frac{(z-z_{0}-\mu Et)^{2}}{2\sigma^{2}}]},$

where $\sum N$ represents the total number of carriers, $\tau$ denotes the
temporal attenuation coefficient of the laser,
$\sigma^{2}=\sigma_{0}^{2}+2Dt$, and $\sigma_{0}$ represents the standard
deviation of the initial Gaussian distribution of carriers. The convolution
asterisk now only affects the time dimension. This analysis is effective for
edge-TCT, as the injected Gaussian beam profile behaves like a Gaussian
function with $\sigma_{0}=w(x)/2$ and $\tau=\Delta t/2\sqrt{2\ln 2}$.

Now consider the elongation of the carrier cluster during drifting. Denote
$D=D(z)$ and $E=E(z)$ as the diffusion constant and the field intensity at the
center of the carrier cluster, respectively. In an infinitesimal time period
$\mathrm{d}t$, the distance between any two carriers would be stretched by the
gradient of the field, which can be represented as $\frac{\mathrm{d}(\mu
E)}{\mathrm{d}z}\mathrm{d}t$. Additionally, the Einstein coefficient $D$ will
change as the carrier cluster center moves, and as a result, the overall
distribution will take this form:

$N(z,t+\mathrm{d}t)=N_{0}\exp{\left\\{-\frac{[z/(1+\frac{\mathrm{d}(\mu
E)}{\mathrm{d}z}\mathrm{d}t)]^{2}}{2(\sigma^{2}+D\mathrm{d}t)}\right\\}}.$
(18)

Therefore,

$\frac{\mathrm{d}\sigma^{2}}{\mathrm{d}t}=\frac{k_{B}T}{q}\mu+2\sigma^{2}\frac{\mathrm{d}(\mu
E)}{\mathrm{d}z}.$ (19)

Thus, we could integrate (the foot marks $i$ and $e$ represent the initial and
end points of the cluster center’s path) and obtain

$\sigma^{2}=\frac{(\mu_{e}E_{e})^{2}}{(\mu_{i}E_{i})^{2}}\left[\int_{z_{i}}^{z_{e}}{\frac{k_{B}T}{q}\frac{\mathrm{d}z}{E}\frac{(\mu_{i}E_{i})^{2}}{(\mu
E)^{2}}+\sigma_{0}^{2}}\right],$ (20)

and if the variation of mobility $\mu$ could be ignored and a uniformly doped
bulk is considered (the electric field intensity of which could be calculated
by $E(z)=E(z_{e})-qN_{eff}(z-z_{e})/\varepsilon$), it is revealed that

$\sigma^{2}=\frac{E_{e}^{2}}{E_{i}^{2}}\left[\sigma_{0}^{2}+\frac{k_{1}}{2}\left(1-\frac{E_{i}^{2}}{E_{e}^{2}}\right)\right],$
(21)

where $k_{1}=(k_{B}T/q)\cdot(\varepsilon/qN_{eff})=2.8\rm\ {\mu m}^{2}$.

Now consider the multiplication process in LGAD. Notice that from then on,
only the carriers leading the avalanche will be considered, as due to the
device design, they often have a larger ionization rate. To simplify
calculations, it could be assumed that during the avalanche, the diffusion
speed of the carrier cluster is slow enough for the spatial distribution to be
treated as Gaussian with a variance that does not change over time. Thus, to
count the carriers that reached the gain layer from a given cluster, the
relation 17 can be deduced into

$N(z_{e},t)=\frac{\sum
N}{\sqrt{2\pi[\tau^{2}+(\sigma/v_{e})^{2}]}}\exp{\left[-\frac{(t-t_{0})^{2}}{2\tau^{2}+2(\sigma/v_{e})^{2}}\right]},$
(22)

where $z_{e}$ refers to the bottom of the gain layer, $v_{e}=\mu E(z_{e})$
represents the drifting velocity of the carriers at that location, and
$t_{0}=\int_{z_{i}}^{z_{e}}dz/v$.

Consider carriers, the number of which is $\mathrm{d}N=N\mathrm{d}t$, passing
through the gain layer within a small time period $\mathrm{d}t$. As a result,
$M\mathrm{d}N$ ionized carrier pairs are created, which then start to drift
and generate a signal. As only carriers with opposite charges can continue to
drift and contribute to the signal while crossing the space charge region, the
amplitude will increase by $M\mathrm{d}Nqv_{e}^{\prime}/d$ (where the prime
denotes carriers with a different charge sign, and $d$ represents the
thickness of the detector). Therefore, the time derivative of the total gain
signal, $\mathrm{d}i_{gain}/\mathrm{d}t$, is proportional to $N$, and the
maximum of the derivative is also proportional to the maximum of $N$. This is
because the gain signal $i_{gain}$ dominates the rising edge of the total
signal $i_{q}$ by

$\frac{\mathrm{d}i_{q}}{\mathrm{d}t}\Big{|}_{max}=\frac{k_{2}\sum
N}{\sqrt{\tau^{2}v_{e}^{2}+\sigma^{2}}},$ (23)

where $k_{2}=(Mq/\sqrt{2\pi})/(v_{e}v^{\prime}_{e}/d)=2.30\times
10^{-7}\rm\frac{A\cdot\mu m}{ns}.$ Notice that for non-irradiated detectors,
$k_{2}$ could be completely determined by the wavelength of the laser and
properties of the detector, since $v_{e}$ is a function of $E_{e}$, and the
value of $E_{e}$ can be estimated by $E_{e}\approx 2(V_{GD}-V_{FD})/d$ where
$V_{GD}$ and $V_{FD}$ represent the full depletion voltages of the gain layer
and the active area.

Considering the distribution of edge-TCT induced carriers described by
Equation 4, in Equation 20 or 21, $\sigma_{0}$ will be replaced by the beam
waist $w(x)/2$. Additionally, it is necessary to integrate along the laser
beam direction $x$ in order to evaluate the change in beam waist. Plus, to
account for attenuation from Equation 5, within the range $(x,x+\mathrm{d}x)$,
the total number of carriers is $(E_{p}/h\nu)\exp{(-\alpha
x)}\alpha\mathrm{d}x$. Thus,

$\frac{\mathrm{d}i_{q}}{\mathrm{d}t}\Big{|}_{max}=\int_{0}^{x_{max}}\frac{k_{2}(E_{p}/h\nu)\exp{(-\alpha
x)}\alpha\mathrm{d}x}{\sqrt{\tau^{2}v_{e}^{2}+\frac{E_{e}^{2}}{E_{i}^{2}}[\frac{w_{0}^{2}(x)}{4}+\frac{k_{1}}{2}(1-\frac{E_{i}^{2}}{E_{e}^{2}})]}}.$
(24)

For facilities utilized in this work, since $\tau v_{e}$ is a relatively large
term (approximately 15 $\mu\rm m$) and $\frac{E(z_{e})}{E(z)}$ is less than
1.1, the result should be approximately constant. Though, we could expect the
diffusion profile method to be more valuable and provide more information in
the case that $v_{e}$ is relatively small.

The simulated and experimental diffusion profiles of LGAD are shown in Figure
7(b), along with the corresponding theoretical values. The reliability of the
diffusion profile method for LGAD has been proven, as it accurately predicts
and observes absolute results.

##  6 Conclusion

The present study has successfully established and tested the diffusion
profile method for measuring the bulk field in LGAD (low-gain avalanche
detector). This methodology aims to use thin detectors with a thickness
approximately ten times the FWHM (full width at half maximum) of the laser
beam. By calculating the temporal derivative of the waveform, valuable
information about the diffusion process can be obtained, as well as its impact
on the flattening of the carrier cluster. The methodology is expected to yield
immediate results by using directly measurable parameters, eliminating the
need for normalization.

This study effectively conducted a simulation of the edge-TCT (edge transient
current technique) on non-irradiated silicon LGAD using the particle detector
simulation package RASER. The simulated waveform was used to fit the data.
Subsequently, an edge-TCT scan was generated by comparing the experimental
results with HPK Type 3.1 detectors (p-i-n and LGAD) to determine its
reliability. The comparison yielded successful results in terms of amplitudes,
rise times, and charge collections, thereby confirming the accuracy of the
simulation in edge-TCT.

A comparable simulation procedure could be implemented to examine the
performance of detectors made from alternative materials, such as silicon
carbide, with a specific focus on their resistance to radiation when exposed
to a particular irradiation flux. Based on prior research [7], simulation can
also be used to investigate TCT tests on detectors with non-planar structures
[32].

## Acknowledgement

This work has been supported by the China Postdoctoral Science Foundation (No.
2022M710085), the State Key Laboratory of Particle Detection and Electronics
(No. SKLPDE-KF-202313), and the Natural Science Foundation of Shandong
Province Youth Fund (No. ZR202111120161). We have completed this work as part
of the CERN RD50 collaboration framework, with the experiment conducted at the
CERN-SSD laboratory. We would like to thank Marcos Fernández García and
Michael Moll for their support during the experiment, as well as the TRACS
team for their valuable cooperation and discussions.

## References

  * [1] X. Shi _et al._ , “raser · PyPI,” https://pypi.org/project/raser/, [accessed June 25, 2023].
  * [2] H. P. K.K., https://www.hamamatsu.com, [accessed June 25, 2023].
  * [3] M. Moll _et al._ , “CERN RD50 - Radiation hard semiconductor devices for very high luminosity colliders,” https://rd50.web.cern.ch, [accessed June 25, 2023].
  * [4] CIVIDEC, “cividec instrumentation - CVD Diamond Technology applications,” https://cividec.at/electronics-C2-HV.html, [accessed June 25, 2023].
  * [5] KEYSIGHT, “DSO9254A Oscilloscope: 2.5 GHz, 4 Analog Channels [Discontinued] — Keysight,” https://www.keysight.com/us/en/product/DSO9254A/oscilloscope-2-5ghz-4-analog-channels.html, [accessed July 16, 2023].
  * [6] T. Yang _et al._ , “Time Resolution of the 4H-SiC PIN Detector,” _Frontiers in Physics_ , vol. 10, 2022. [Online]. Available: https://www.frontiersin.org/articles/10.3389/fphy.2022.718071
  * [7] Y. Tan _et al._ , “Timing Performance Simulation for 3D 4H-SiC Detector,” _Micromachines_ , vol. 13, no. 1, 2022. [Online]. Available: https://www.mdpi.com/2072-666X/13/1/46
  * [8] R. Van Overstraeten and H. De Man, “Measurement of the ionization rates in diffused silicon p-n junctions,” _Solid-State Electronics_ , vol. 13, no. 5, pp. 583–608, 1970. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0038110170901395
  * [9] S. Reggiani, M. Valdinoci, L. Colalongo, M. Rudan, G. Baccarani, A. D. Stricker, F. Illien, N. Felber, W. Fichtner, and L. Zullino, “Electron and hole mobility in silicon at large operating temperatures. i. bulk mobility,” _IEEE Transactions on Electron devices_ , vol. 49, no. 3, pp. 490–499, 2002.
  * [10] G. Kramberger _et al._ , “Investigation of Irradiated Silicon Detectors by Edge-TCT,” _IEEE Transactions on Nuclear Science_ , vol. 57, no. 4, pp. 2294–2302, 2010.
  * [11] ——, “Electric ﬁeld modeling in heavily irradiated silicon detectors based on Edge-TCT measurements. ”Why p-type is Better than n-type” or Electric Field in Heavily Irradiated Silicon Detectors,” _PoS_ , vol. Vertex2012, p. 022, 2013. [Online]. Available: https://cds.cern.ch/record/2634361
  * [12] B. Schmidt, “The high-luminosity upgrade of the LHC: Physics and technology challenges for the accelerator and the experiments,” _Journal of Physics: Conference Series_ , vol. 706, no. 2, p. 022002, apr 2016. [Online]. Available: https://dx.doi.org/10.1088/1742-6596/706/2/022002
  * [13] V. Izzo, “ATLAS Upgrades,” CERN, Geneva, Tech. Rep., 2021. [Online]. Available: https://cds.cern.ch/record/2732959
  * [14] R. Longo, “Joint ATLAS/CMS ZDC upgrade project for the High Luminosity LHC,” CERN, Geneva, Tech. Rep., 2022. [Online]. Available: https://cds.cern.ch/record/2825363
  * [15] L. Soffi, “Physics opportunities with a MIP Timing Detector in CMS at HL-LHC,” CERN, Geneva, Tech. Rep., 2022. [Online]. Available: https://cds.cern.ch/record/2840966
  * [16] G. Calderini, “The ATLAS ITk detector for High Luminosity LHC upgrade,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 1040, p. 167048, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900222004715
  * [17] H. F. W. Sadrozinski _et al._ , “Sensors for ultra-fast silicon detectors,” _Nucl. Instrum. Methods A_ , vol. 765, pp. 7–11, 2014. [Online]. Available: http://dx.doi.org/10.1016/j.nima.2014.05.006
  * [18] Z. Galloway _et al._ , “Properties of HPK UFSD after neutron irradiation up to 6e15 n/cm2,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 940, pp. 19–29, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900219306278
  * [19] X. Yang _et al._ , “Layout and performance of HPK prototype LGAD sensors for the High-Granularity Timing Detector,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 980, p. 164379, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900220307762
  * [20] G. Giacomini _et al._ , “Fabrication of Different LGAD-Based Devices at BNL,” in _2021 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_ , 2021, pp. 1–4.
  * [21] M. O. Wiehe, “Development of a two-photon absorption - tct system and study of radiation damage in silicon detectors,” Ph.D. dissertation, Freiburg U., Freiburg U., 2021.
  * [22] B. J. Baliga, _Fundamentals of Power Semiconductor Devices_. Springer, 2019.
  * [23] J. Calvo _et al._ , “TRACS: A multi-thread transient current simulator for micro strips and pad detectors,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 917, pp. 77–85, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900218317856
  * [24] K. Seeger, _Semiconductor physics: An introduction_. Springer, 2011.
  * [25] K. Thyagarajan and A. Ghatak, _Lasers: Fundamentals and Applications_ , ser. Graduate Texts in Physics. Springer US, 2010. [Online]. Available: https://books.google.com.hk/books?id=z13wEOBwn1wC
  * [26] S. Spannagel _et al._ , “Allpix Squared - Generic Pixel Detector Simulation Framework,” https://doi.org/10.5281/zenodo.6387859, [accessed June 25, 2023].
  * [27] F. Cenna _et al._ , “Weightfield2: A fast simulator for silicon and diamond solid state detector,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 796, pp. 149–153, 2015, proceedings of the 10th International Conference on Radiation Effects on Semiconductor Materials Detectors and Devices. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900215004842
  * [28] M. S. Alnaes _et al._ , “The FEniCS project version 1.5,” _Archive of Numerical Software_ , vol. 3, 2015.
  * [29] V. Eremin _et al._ , “Development of transient current and charge techniques for the measurement of effective net concentration of ionized charges (Neff) in the space charge region of p-n junction detectors,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 372, no. 3, pp. 388–398, 1996. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0168900295012958
  * [30] N. Moffat _et al._ , “Low Gain Avalanche Detectors (LGAD) for particle physics and synchrotron applications,” _Journal of Instrumentation_ , vol. 13, no. 03, p. C03014, mar 2018. [Online]. Available: https://dx.doi.org/10.1088/1748-0221/13/03/C03014
  * [31] V. Eremin _et al._ , “The origin of double peak electric field distribution in heavily irradiated silicon detectors,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 476, no. 3, pp. 556–564, 2002, proc. of the 3rd Int. Conf. on Radiation Effects on Semiconductor Materials, Detectors and Devices. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900201016424
  * [32] L. Diehl _et al._ , “Comparison of the time resolution of unirradiated and irradiated LGADs and 3D sensors,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 1046, p. 167691, 2023. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168900222009834
  * [33] S. Selberherr _et al._ , “The evolution of the MINIMOS mobility model,” _Solid-State Electronics_ , vol. 33, no. 11, pp. 1425–1436, 1990. [Online]. Available: https://www.sciencedirect.com/science/article/pii/003811019090117W


# Modeling Current and Future High-Cadence Surveys of Repeating FRB
Populations

Kyle McGregor Wesleyan University Department of Astronomy and Van Vleck
Observatory, 96 Foss Hill Dr., Middletown CT 06459, USA  Department of Physics
and Astronomy, West Virginia University, P.O. Box 6315, Morgantown, WV 26506,
USA  Duncan R. Lorimer Department of Physics and Astronomy, West Virginia
University, P.O. Box 6315, Morgantown, WV 26506, USA  Center for Gravitational
Waves and Cosmology, West Virginia University, Chestnut Ridge Research
Building, Morgantown, WV 26505, USA

###### Abstract

In recent years, the CHIME (Canadian Hydrogen Intensity Mapping Experiment)
interferometer has revealed a large number of Fast Radio Bursts (FRBs),
including a sizable population that demonstrates repeating behavior. This
transit facility, employing a real-time FRB search pipeline, continually scans
the sky with declinations between $-10^{\circ}$ and $90^{\circ}$ for events
with fluences $\gtrapprox 0.4$ Jy ms. We simulate a population of repeating
FRBs by performing Monte Carlo simulations of underlying source populations
processed through a mock CHIME/FRB observing pipeline. Assuming intrinsic
repeater rates follow a Poisson distribution, we test assumptions about the
burst populations of the repeater sample, and construct models of the FRB
sample assuming various cosmological distributions. We infer the completeness
of CHIME/FRB observations as a function of observing cadence and redshifts out
to 0.5. We find that, if all simulated bursts have a fixed Poisson probability
of repetition over their integrated time of observation, repeating burst
detections across comoving volume should continue to grow near linearly on the
order of decades. We predict that around 170 of the current CHIME/FRB one-off
sources will ultimately repeat. We also make projections for FRB repeaters by
future facilities and demonstrate that the number of repeaters they find could
saturate on a $\sim 3$ yr timescale.

Unified Astronomy Thesaurus concepts: Radio transient sources (2008)

††software: astropy (Astropy Collaboration et al., 2013),  scipy (Virtanen et
al., 2020),  pygedm (Price et al., 2021)

##  1 Introduction

Fast Radio Bursts (FRBs) are among the most enigmatic extragalactic sources of
radiation. First discovered in 2007 using archival data from the Parkes
(Murriyang) telescope (Lorimer et al., 2007) these dispersed highly-energetic
millisecond-scale bursts of radio emission have been shown to be consistent
with a cosmological origin (Thornton et al., 2013). While no conclusive
progenitor mechanism for FRB emission has been widely accepted, modeling shows
that directional magnetar flares have energy budgets similar to observed burst
energies (see, e.g., Popov et al., 2018).

With 672 one-off bursts currently known111For a list, see
https://blinkverse.alkaidos.cn. (Xu et al., 2023), a recent data release from
the Canadian Hydrogen Intensity Mapping Experiment (CHIME) collaboration has
more than doubled the number of known repeating FRBs, currently standing at 63
(CHIME/FRB Collaboration et al., 2023). Searches led by the CHIME/FRB
collaboration using this facility are carrying out a census of FRBs. As a
transit instrument that each day observes the 200 square degrees flanking the
north-south meridian of the sky at declinations between –10 and 90 degrees,
CHIME repeatedly scans the sky and is therefore well suited to probe the
repeating FRB population.

With a small known sample size compared to the inferred daily incidence of
bursts from across the universe (expected to be on the order of thousands per
day; Thornton et al., 2013), a Monte Carlo simulation is a powerful tool to
infer properties of the underlying FRB population(s). It is also important to
determine the limits of our facilities and find how many bursts can be found
in a given cadence, as well as account for survey selection effects that may
influence completeness. We present a model for an observed repeater population
constructed from first principles, which can provide a prediction for future
detection rates by CHIME/FRB. From this, we can determine a model census of
bursts whose properties can be compared to the known sample revealed in
(CHIME/FRB Collaboration et al., 2023).

Inspired by the latest release of the repeater sample (CHIME/FRB Collaboration
et al., 2023), we apply a Monte Carlo approach in this paper. In Section 2 we
discuss how we build a model source population, in Section 3 we discuss our
procedure for generating a mock CHIME/FRB observing campaign on this
population. We compare these modeled populations to the observed population in
Section 4, and discuss the relevance to previous and future work in Section 5
and present our conclusions in Section 6.

##  2 Model Source Properties

In constructing a synthetic population of repeaters, we must first consider
properties intrinsic to each source. These parameters include right ascension
and declination on the sky ($\alpha$, $\delta$), redshift ($z$), mean
luminosity ($L_{0}$), intrinsic pulse width ($w_{\rm int}$), and host-galaxy
dispersion measure ($\rm{DM}_{Host}$).

As FRBs are a cosmological population, we assume a uniform distribution of
locations on the celestial sphere. We draw from a uniform random deviate
between 0–24 hr in right ascension. In declination space, however, following
Chawla et al. (2022), we draw uniformly across
$-10^{\circ}<\delta<90^{\circ}$. The cosine dependence in solid angle for an
isotropic population cancels the inverse cosine dependence in exposure. While
not representing a uniform distribution on the celestial sphere, this does
give a distribution proportional to the product of the declination-dependent
sky coverage and exposure for CHIME/FRB which more accurately represents the
underlying position distribution of detected sources (for details, see Chawla
et al., 2022).

It is currently unclear what the underlying redshift distribution of FRBs is
in general, and there is debate as to the distribution probed by CHIME/FRB
(see, e.g., Zhang & Zhang, 2022; Shin et al., 2023). In addition, the DM
distributions for CHIME/FRB repeaters appears to be significantly different to
those of the one-off bursts (CHIME/FRB Collaboration et al., 2023), which
leads to further uncertainty in their underlying redshift distributions. For
our models of repeating FRBs, we explore two scenarios for the redshift
distributions: a population following a constant density in comoving volume
($V_{c}$) and a population following the cosmic star formation rate (SFR;
Madau & Dickinson, 2014). For the former case, the probability density in
redshift

$P_{\rm comoving}(z)\propto\left(\frac{1}{1+z}\right)\frac{dV_{c}}{d\Omega
dz},$ (1)

where $\Omega$ is the solid angle and $\frac{dV_{c}}{d\Omega dz}$ is the
differential comoving volume element. For the latter case,

$P_{\rm
SFR}(z)\propto\left(\frac{(1+z)^{2.7}}{1+\left[\frac{(1+z)}{2.9}\right]^{5.6}}\right)\,\left(\frac{1}{1+z}\right)\,\frac{dV_{c}}{d\Omega
dz}\cdot$ (2)

In both cases, as we discuss later, we restrict our maximum redshifts probed
by the simulations to be $z_{\rm max}=0.5$. More distant sources are not
required to account for the currently observed sample of repeaters.

For the mean luminosity of each burst, $L_{0}$, we draw from a Schechter
function employing a fit from Shin et al. (2023). We impose bounds between
$10^{39}$ and $10^{45}$ erg s-1. For a 1 ms pulse, this is consistent with the
inferred luminosity limits of observed bursts, consistent with the burst
energy limits used by Chawla (2022). The Schechter luminosity distribution is
commonly used to describe the energetics of extragalactic sources and takes
the form

$P_{L_{0}}(L_{0})dL_{0}=\left(\frac{L_{0}}{L_{*}}\right)^{\alpha+1}\exp{\left(\frac{L_{0}}{L_{*}}\right)}d(\log
L_{0}),$ (3)

where the characteristic luminosity $L^{*}=2\times 10^{44}$ erg s-1 and power
law slope $\alpha=-1.3$. Given that the current CHIME/FRB repeaters are known
to vary in intensity (CHIME/FRB Collaboration et al., 2023), we dither the
daily burst luminosities using a scheme described in Section 3.

In simulating intrinsic pulse widths, we draw from a lognormal distribution as
given in Luo et al. (2020). The dispersion measure (DM) of a burst is the
integral of electron number density over the line of sight to the burst, which
will include contributions from each source of propagation en route to Earth
and thus serves as an observable proxy for distance. For our simulations,
following previous authors (e.g., Luo et al., 2018), we computed the observed
DM,

$\rm
DM_{obs}=DM_{MW}(\textit{l},\textit{b})+DM_{IGM}(\textit{z})+\frac{DM_{host}}{1+\textit{z}},$
(4)

with the Galactic component, DMMW, from the YMW16 model of Milky Way electron
density (Yao et al., 2017), the intergalactic medium (IGM) component, DMIGM,
from the Macquart ($z$-DM) relation (James et al., 2022), and the repeater
host galaxy contribution, DMhost, from a lognormal random deviate parametrized
by Mo et al. (2023).

We draw a catalog of $10^{5}$ mock repeaters using Monte Carlo methods in
scipy. This gives us a base of source repeater progenitors over which we can
iterate to simulate an observation campaign with CHIME, the procedure for
which is described in the next section.

##  3 Model CHIME Observing Campaign

While CHIME is a complex instrument, responsive but highly sensitive to on-
site conditions that are prohibitive to model efficiently (see, e.g., Andersen
et al., 2023), we can infer the observability of a given event by imposing a
fluence threshold of 0.4 Jy ms. This lower limit is consistent across the
CHIME Catalog 1 data (CHIME/FRB Collaboration et al., 2021) and the new
repeater data (CHIME/FRB Collaboration et al., 2023), and used in similar
population syntheses (Chawla, 2022).

Figure 1: Fluence (top) and flux (bottom) distributions for the CHIME/FRB
repeater sample (green histograms) showing the drop in the number of bursts
around the threshold. Dashed and solid lines show polynomial fits to these
tails in log and linear spaces, respectively. Normalization of these fits on
the right to a probability scale allows us to implement these sensitivity roll
offs in our simulated FRBs (see text). Figure 2: Observed CHIME/FRB fluence
distributions compared to those from our model population. The simulated
distribution (orange) shows the constant start formation model with the simple
0.4 Jy ms fluence threshold applied. The selection-corrected distribution
(blue) shows just those FRBs which survived the Bernoulli trials in flux and
fluence (see text).

Our mock observation pipeline first iterates over each of the $10^{5}$
simulated repeater and determines the daily burst incidence using a Poisson
random deviate. Due to the shape of CHIME’s primary beam, higher declinations
will have a greater exposure over CHIME’s operations, resulting in a more
complete sample at high declinations. This declination dependence has been
modeled in similar studies (see James, 2023). We choose to weigh the daily
Poisson burst incidence by this exposure in our forward modeling, such that
higher declinations are more likely to sound off on a given day to correct for
the declination bias. The mean of the probability mass function for the number
of observed bursts each day over the four CHIME cylinders,

$\lambda=\frac{4\times\textrm{FWHM}_{E-W}}{360^{\circ}\cos(\delta)}\,\,\Gamma(R_{\rm
min},R_{\rm max},\gamma),$ (5)

with $\textnormal{FWHM}_{\rm E-W}=0.32^{\circ}$ the full width at half maximum
of each of synthesized beam (Chawla, 2022). This corresponds to the product
$\lambda=(\textrm{time in beam})\times(\textrm{bursts day}^{-1})$ for each day
of the simulation. Following James (2023), $\Gamma(R_{\rm min},R_{\rm
max},\gamma)$ is a power-law random deviate ranging between $R_{\rm
min}=10^{-3}$ and $R_{\rm max}=10^{-1}$ $\rm bursts~{}day^{-1}$ with a power-
law index $\gamma=-2.2$. We choose $R_{\rm min}$ to be approximately the
reciprocal of the current survey length (3.3 yr) and models with $R_{\rm
max}>0.1$ day-1 produce too many observable events. These are defined in terms
of the total survey cadence, rather than the time in beam. In addition, we do
not assume any correlation between burst rate and intrinsic luminosity,
$L_{0}$, over the range of $L_{0}$ we model.

The observed pulse width, $w_{\rm obs}$, includes contributions from sources
of propagation and instrumental parameters. It is given by the quadrature
addition

$w_{\rm obs}=\sqrt{[(1+z)w_{\rm int}]^{2}+w_{\rm scat}^{2}+w_{\rm
samp}^{2}+w_{\rm DM}^{2}}\,,$ (6)

where $w_{\rm int}$ is the intrinsic pulse width and the $(1+z)$ factor
accounts for time dilation, $w_{\rm scat}$ is the scatter broadening, $w_{\rm
samp}$ is the data sampling interval and $w_{\rm DM}$ is the dispersion
broadening across a finite frequency channel. To model the intrinsic width,
following Luo et al. (2020), we draw from a log normal distribution so that
log${}_{10}\,(w_{\rm int}/{\rm ms})=$ normal($\mu_{W},\sigma_{W}$). Here the
normal distribution has a mean $\mu_{W}=0.2$ and standard deviation
$\sigma_{W}=0.33$. The CHIME/FRB sampling time $w_{\rm samp}=0.98$ ms (CHIME
Collaboration et al., 2022). The scattering timescale is calculated alongside
the $\rm DM_{MW}$ and $\rm DM_{IGM}$ contributions with the pygedm package
(Price et al., 2021) using the electron density model developed by Yao et al.
(2017).

As mentioned above, some variation in the flux density of repeaters across
observations is noted in the CHIME sample, which implies a mechanism for
differing luminosity over consecutive bursts from the same repeater. While the
sample size for most observed bursts is much too small to easily parametrize,
we assume a normal distribution with a standard deviation $\sigma=0.1L_{0}$.
Using Eq. 8 of Macquart & Ekers (2018), we find that the observed fluence

$F_{\nu}=\frac{(1+z)\,w_{\rm int}}{4\pi D_{L}^{2}\Delta\nu}\,\,{\rm
normal}(L_{0},0.1L_{0}),$ (7)

where $D_{L}$ is the luminosity distance and the bandwidth for CHIME/FRB,
$\Delta\nu=400$ MHz. Notably here, for simplicity, we are being agnostic about
any dependence of the fluence with frequency and implicitly assume that the
FRBs are flat spectrum sources over the frequency ranges explored here (i.e.,
400 MHz out to $(1+z_{\rm max})\,800$ MHz = 1200 MHz).

As noted by Cordes & McLaughlin (2003), the signal-to-noise threshold for FRB
detection can be computed through the radiometer noise considerations.
Originally, we followed this approach and implemented a discretization of
CHIME/FRB’s beam to infer the observability of a given FRB using the
radiometer equation. This method is unphysical in that it assumes a fiducial
nature to CHIME’s system-equivalent flux density, which in reality is
calibrated in real-time as part of CHIME/FRB’s bonsai FRB search algorithm.
Following the detailed sensitivity determinations by Chawla (2022), however,
after iterating over each day for each repeater in the population, we accept a
burst as observable if its fluence is above the 0.4 Jy ms threshold. This
gives us a population theoretically observable by CHIME, though selection
corrected for various incompleteness factors compared to the CHIME survey
populations.

To accurately model CHIME’s observing campaign, we must reintroduce these
incompleteness factors which become important close to the detection
threshold. It is known that CHIME is observationally biased against low-
fluence events, which can only be detected at the highest gain regions of each
synthesized beam, whereas higher fluence events can be readily detected across
the primary beam and sidelobes (Lin et al., 2023). We account for this bias by
weighting the probability of observation for each low-fluence and low-flux
event with an independent Bernoulli trial, which will return a “success” or
“failure” based on the exhausted probability of detection at each fluence and
flux in the over-represented regime. We employ a polynomial fit for the drop-
off following the empirical cutoff curve in the observed fluence/flux
distributions between the 0.4 Jy ms lower fluence limit and the 5 Jy ms $95\%$
completeness threshold identified in CHIME/FRB Collaboration et al. (2021) .
This is done during the filtering step after the mock observed catalog is
collected in order to aid computational efficiency.

The fourth-degree polynomial fits to the left edges of the fluence and flux
distributions are shown in Fig. 1. We report the normalized fits as
$p(F_{\nu})=-0.002(F_{\nu}-5)^{4}+1$ for fluence and
$p(S_{\nu})=-240.8S_{\nu}^{4}+174.4S_{\nu}^{3}-30.85S_{\nu}^{2}+2.46S_{\nu}-0.045$
for flux density where the units of fluence and flux density at Jy ms and Jy,
respectively. In Fig. 2 we show the application of this approach to the
constant star formation model. The better agreement between the
“Observability” and “Completeness” thresholds is clearly evident.

##  4 Results

We ran the simulations above for a total of 7300 days (20 years), extracting
the sample at 1095 days (three years) to match the timescale of the CHIME
sample (CHIME/FRB Collaboration et al., 2023). Fig. 3 show cumulative
distributions for the observed and model samples for the constant star
formation and constant volume redshift distributions, respectively. For each
model, we show the fluence, pulse width, DM, and declination distributions.
While we do not take the approach of fitting models to the data in this paper,
it is useful to quote quantitative metrics of comparison between the model and
observed samples. Table 1 shows the results of the two-sided Kolmogorov-
Smirnov (KS) tests (Kolmogorov, 1933; Smirnov, 1948) we ran on these
distributions.

Table 1: KS-test results for the two redshift models considered in this study. For each model, we provide the KS statistic and associated $p$-value for fluence, observed pulse width, DM and declination, respectively.  Model |  | $F_{\nu}$ | $w_{obs}$ | DM | Decl.
---|---|---|---|---|---
Eq. 1 | KS statistic | $0.383$ | $0.189$ | $0.112$ | $0.154$
|  $p$-value | $2\times 10^{-8}$ | $0.030$ | $0.398$ | $0.112$
Eq. 2 | KS statistic | $0.406$ | $0.185$ | $0.104$ | $0.146$
|  $p$-value | $3\times 10^{-9}$ | $0.035$ | $0.565$ | $0.147$

Figure 3: Cumulative distributions for the redshift distribution tracing the
constant volume redshift distribution (Eq. 1; upper row) and the cosmic star
formation rate (Eq. 2; lower row). Green lines show the observed distributions
from CHIME/FRBs. The blue lines are our simulated observable populations.

##  5 Discussion

The goal of this work is to constrain some of the properties of the FRB
repeater population as seen through regular observations with the CHIME/FRB
survey described by CHIME/FRB Collaboration et al. (2023). In doing this, we
have tried to strike a balance between a simple analytic model and a very
detailed simulation which would also include Markov Chain Monte Carlo
techniques for parameter estimation. In many respects, we have taken an
approach similar to earlier studies of the pulsar population (see, e.g.,
Faucher-Giguère & Kaspi, 2006) where plausible models are developed that can
mimic the main results. In the sections below, we discuss the main findings
from this approach.

###  5.1 Which model is better?

As can be seen in Fig. 3, and summarized quantitatively in Table 1, the
agreement between our model and observed FRB populations is generally good and
we briefly comment on each of the distributions considered. While the model
with the redshift distribution tracing the cosmic star formation matches the
DM distribution slightly better than the constant volume model, the KS tests
indicate that the two models are statistically very similar in their match to
the observed CHIME/FRB repeater DM distribution. This is also reflected in the
other distributions and KS test scores and is generally in lines with our
expectations, since the redshift distributions we considered have a similar
behaviour in the range $0<z<0.5$. For the fluence distribution, we note that
our model fluences are systematically higher than those observed by CHIME/FRB.
We did not seek to optimize this further, given that the fluences in the
CHIME/FRB catalog are lower limits due to uncertainties in the true source
position within the CHIME beam. In the remainder of the discussion, we use the
cosmic star formation redshift distribution (Eq. 2) as our reference model.

###  5.2 Repeater detections versus time

We have focused our attention on the cumulative distributions in Fig. 3, but
note that our simulations naturally predict the detections as a function of
time. As an example, to compare with Fig. 3 of CHIME/FRB Collaboration et al.
(2023), in Fig. 4 we show a set of randomly selected repeaters from our
reference model as a function of declination versus time. While we see a
qualitative similarity to Fig. 3 of CHIME/FRB Collaboration et al. (2023), we
note that (as anticipated from Fig. 3), our modeling approach results in an
overabundance of sources at higher declinations than actually observed. As
seen by CHIME/FRB Collaboration et al. (2023), it is notable that the higher
declination sources benefit from a longer exposure time due to the unique
optics of CHIME resulting in more bursts observed. Further modeling, including
different burst rate distributions would be a logical extension of this work.

Figure 4: Declination versus time for 163 selected repeaters from the
simulated population. Cyan represents bursts with $F_{\nu}<4$ Jy ms, magenta
$4\textrm{ Jy ms}<F_{\nu}<40$ Jy ms, and yellow $F_{\nu}>40$ Jy ms.

One aspect we do consider in this paper, however, are some straightforward
predictions for the number of repeaters we expect with future CHIME/FRB
observations. The power-law burst rate distribution we have considered from
James (2023) produces a good match for the data through CHIME/FRB
Collaboration et al. (2023), which to date has reported a near-linear pattern
of detection incidence. After scaling the number of sources to match the
number of total number of cumulative detections by CHIME/FRB to date, our
model predicts this almost linear trend will continue for the next decade, in
agreement with James (2023). Any drop in this overall detection rate within
the next few years would clearly challenge this model.

Figure 5: Cumulative repeater incidence for CHIME/FRB along with our model
projection for the next two decades. We scale the modeled detection incidence
to the first 3.3 years of CHIME/FRB results. Inset: zoomed-in view of the
current sample (CHIME/FRB Collaboration et al., 2023) which shows an
approximately linear trend.

###  5.3 One-off FRBs in the observed population

An interesting question that our modeling approach can address is what the
number of apparent one-off sources in the current CHIME sample are actually
repeaters. Since our modeling is exclusively assuming a repeating population,
we automatically keep track of any simulated FRBs for which only one pulse is
recorded over any given timescale. For the constant star formation case, we
find that, over a three-year period, for every repeater observed, there are
2.8 apparently one-off sources. In other words, for the 60 repeaters currently
observed by CHIME/FRB Collaboration et al. (2023), our model predicts that
around 170 of the currently one-off sources in the CHIME/FRB catalog will
repeat. Based on the projection shown in Fig. 5, we anticipate that these
sources are very likely to repeat within 5 yr.

###  5.4 Predictions for other instruments

Our simulations can also be used to provide interesting predictions for
hypothetical FRB search campaigns with two forthcoming next-generation radio
observatories, the Deep Synoptic Array (DSA-2000) and the Canadian Hydrogen
Observatory and Radio-transient Detector (CHORD). These facilities, which are
planned for construction and first light within the decade, are designed to
have greater sensitivity than CHIME. DSA-2000 is projected to have a field of
view FoV = 10.6 deg2, with a fluence detection threshold for a 1 ms burst of
0.03 Jy ms and a bandwidth of 850 MHz in the 1–2 GHz band (Hallinan et al.,
2019). CHORD is slated to be the direct successor to CHIME, with FoV = 130
deg2, a fluence detection threshold for a 1 ms burst of 0.1 Jy ms, and a
bandwidth of 1200 MHz in the 0.3–1.4 GHz band (Vanderlinde et al., 2019).

Figure 6: Cumulative detection predictions for DSA2000 and CHORD compared to
CHIME/FRB prediction from Fig. 5. For display purposes, the model projections
have been smoothed by fitting their cumulative incidences to gamma
distributions. The upper and lower limits are chosen to reflect two extremes
of observational strategies for observing a specific field of view.

Using these specifications, to illustrate the impacts of these sensitivity
improvements, we proceed with two cases acting as upper and lower limits for
the daily time in beam (exposure) for each source for each instrument,
respectively. As an upper limit, we consider a sample of sources constantly
observed over the entire simulated campaign, corresponding to a circumpolar or
never-set population always within the array’s field of view, giving a daily
exposure of 1 day. For a lower limit, we assume a population of same sources
at the celestial equator, such that their daily exposure is $\sqrt{\rm
FoV}/360^{\circ}$ day. To enforce the completeness dropoff approaching the
fluence limit, we scale the fit for $p(F_{\nu})$ assuming an order of
magnitude drop-off between the lower fluence limit and completeness limit
($F_{\rm lower}\sim 10F_{\rm complete}$) for each instrument. This is likely
to be a conservative estimate as these facilities will be more sensitive to
bursts within $z=0.5$ than CHIME.

We show the range between these limiting cases for exposure in Fig. 6. While
we do not have a baseline for the size of population observed at higher
sensitivities, as the scaling of the population discussed in section 5.2 is
unique to a population observed with CHIME’s sensitivity and unique
declination-dependence not shared by these next-generation facilities, we can
still make predictions for the detection rates by these observatories relative
to the simulated population size. Fig. 6 shows the cumulative repeater
detections as a fraction of the total population. Unlike the case for
CHIME/FRB, the more sensitive instruments appear to show a saturation in the
number of repeating FRBs observed starting at around 1000 days.

While the details of the putative surveys with DSA-2000 and CHORD are overly
simplistic, and are restricted to the repeating FRB population with $z<0.5$,
the results shown here highlight the potential for these and other future
facilities with similar sensitivity to probe this population. In the meantime,
we note that our predictions for CHIME are testable in the coming years and
that CHIME’s repeater sample is already a rich resource for studying the FRB
population.

##  6 Conclusions

We have described a simple simulation of the repeating FRB population that is
based on the sample reported recently by CHIME/FRB Collaboration et al.
(2023). We find that this sample of repeaters can be well described by simple
model in which the sources are restricted to redshifts $z<0.5$. Within this
range, we cannot distinguish between a redshift distribution for the
progenitors of FRB repeaters that follows the cosmic star formation, or is
constant in comoving volume. As future CHIME/FRB observations are collected,
as pointed out by James (2023), we anticipate future studies being more
sensitive to the population of FRBs with $z>0.5$.

In spite of the daily cadence of the CHIME/FRB observations, our models
predict that the number of observed repeaters will not saturate significantly
in the coming years and that the sample will grow in size approximately at the
current rate of 20 new repeaters per year. A caveat to this prediction is its
dependence on our assumptions. Future data releases from CHIME/FRB will allow
us to constrain the model parameters/assumptions from this study. Among the
future discoveries our models also predict in the next 5 years are around 170
sources that are currently in the CHIME/FRB sample as “one-off” FRBs.

Our simulation approach has tried to strike a balance between simplicity and
rigor. We have attempted to incorporate the most important aspects of the
population and detection process into our work, but have not tried to fine
tune or do parameter estimation of the model parameters. Further modeling of
the sample which explores the luminosity function, burst rate distribution and
sky exposure are certainly warranted but beyond the scope of the current work.
In particular, we have not accounted for the spectral behavior of the model
FRBs, and the possibility of non-Poissonian burst distributions. Studies of
this nature would be extremely valuable to better understand the growing
population of repeating FRBs.

This research was carried out during a 10-week Research Experience for
Undergraduates (REU) program at West Virginia University. We gratefully
acknowledge the National Science Foundation’s support of the REU under award
number 1950617, as well as the NASA West Virginia Space Grant Consortium for
their support. We thank Emmanuel Fonseca, Clancy James, Marcus Merryfield and
Vicky Kaspi for useful discussions.

## References

  * Andersen et al. (2023) Andersen B. C., et al., 2023, arXiv e-prints, p. arXiv:2305.11302
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * CHIME Collaboration et al. (2022) CHIME Collaboration et al., 2022, ApJS, 261, 29
  * CHIME/FRB Collaboration et al. (2021) CHIME/FRB Collaboration et al., 2021, ApJS, 257, 59
  * CHIME/FRB Collaboration et al. (2023) CHIME/FRB Collaboration et al., 2023, ApJ, 947, 83
  * Chawla (2022) Chawla P., 2022, PhD thesis, McGill University, Canada
  * Chawla et al. (2022) Chawla P., et al., 2022, ApJ, 927, 35
  * Cordes & McLaughlin (2003) Cordes J. M., McLaughlin M. A., 2003, ApJ, 596, 1142
  * Faucher-Giguère & Kaspi (2006) Faucher-Giguère C.-A., Kaspi V. M., 2006, ApJ, 643, 332
  * Hallinan et al. (2019) Hallinan G., et al., 2019, in Bulletin of the American Astronomical Society. p. 255 (arXiv:1907.07648), doi:10.48550/arXiv.1907.07648
  * James (2023) James C. W., 2023, arXiv e-prints, p. arXiv:2306.17403
  * James et al. (2022) James C. W., Prochaska J. X., Macquart J. P., North-Hickey F. O., Bannister K. W., Dunning A., 2022, MNRAS, 509, 4775
  * Kolmogorov (1933) Kolmogorov A. N., 1933, Giornale dell’Istituto Italiano degli Attuari, 4, 83
  * Lin et al. (2023) Lin H.-H., et al., 2023, arXiv e-prints, p. arXiv:2307.05261
  * Lorimer et al. (2007) Lorimer D. R., Bailes M., McLaughlin M. A., Narkevic D. J., Crawford F., 2007, Science, 318, 777
  * Luo et al. (2018) Luo R., Lee K., Lorimer D. R., Zhang B., 2018, MNRAS, 481, 2320
  * Luo et al. (2020) Luo R., Men Y., Lee K., Wang W., Lorimer D. R., Zhang B., 2020, MNRAS, 494, 665
  * Macquart & Ekers (2018) Macquart J. P., Ekers R., 2018, MNRAS, 480, 4211
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, ARA&A, 52, 415
  * Mo et al. (2023) Mo J.-F., Zhu W., Wang Y., Tang L., Feng L.-L., 2023, MNRAS, 518, 539
  * Popov et al. (2018) Popov S. B., Postnov K. A., Pshirkov M. S., 2018, Physics Uspekhi, 61, 965
  * Price et al. (2021) Price D. C., Flynn C., Deller A., 2021, PASA, 38, e038
  * Shin et al. (2023) Shin K., et al., 2023, ApJ, 944, 105
  * Smirnov (1948) Smirnov N., 1948, The Annals of Mathematical Statistics, 19, 279
  * Thornton et al. (2013) Thornton D., et al., 2013, Science, 341, 53
  * Vanderlinde et al. (2019) Vanderlinde K., et al., 2019, in Canadian Long Range Plan for Astronomy and Astrophysics White Papers. p. 28 (arXiv:1911.01777), doi:10.5281/zenodo.3765414
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Xu et al. (2023) Xu J., et al., 2023, Universe, 9, 330
  * Yao et al. (2017) Yao J. M., Manchester R. N., Wang N., 2017, ApJ, 835, 29
  * Zhang & Zhang (2022) Zhang R. C., Zhang B., 2022, ApJ, 924, L14