Geant4-based simulations of charge collection in CMOS Active Pixel Sensors

Geant4 is an object-oriented toolkit for the simulation of the interaction of particles and radiation with matter. It provides a snapshot of the state of a simulated particle in time, as it travels through a specified geometry. One important area of application is the modelling of radiation detector systems. Here, we extend the abilities of such modelling to include charge transport and sharing in pixelated CMOS Active Pixel Sensors (APSs); though similar effects occur in other pixel detectors. The CMOS APSs discussed were developed in the framework of the PRaVDA consortium to assist the design of custom sensors to be used in an energy-range detector for proton Computed Tomography (pCT). The development of ad-hoc classes, providing a charge transport model for a CMOS APS and its integration into the standard Geant4 toolkit, is described. The proposed charge transport model includes, charge generation, diffusion, collection, and sharing across adjacent pixels, as well as the full electronic chain for a CMOS APS. The proposed model is validated against experimental data acquired with protons in an energy range relevant for pCT.


Introduction
Geant4 [1] is an object-oriented toolkit for the simulation of the interactions of particle and radiation with matter. It provides advanced functionality for all the parameters related to detector simulations: geometry and material specifications, definition of particles, physical processes, tracking, event and run management, user interface and visualisation. Although originally developed for High Energy Physics experiments at CERN, it now finds application in many domains of experimental physics, including astrophysics, medical physics or space engineering [2], as particle interactions over a wide range of energies (from 250 eV up to TeV) can be simulated. Geant4 provides a snapshot of the state of a particle at each time instant, while travelling through a specified geometry, and, particularly, permits the energy deposition along the particle path to be determined. However, in most practical cases, simulations are performed to predict the response of a detector in a particular experiments. For this task to be accomplished, simulations need to include charge generation, diffusion, collection and an accurate description of the detector electronics under investigation.
In this paper we report on Geant4-based simulations of charge collection in CMOS Active Pixel Sensors (APSs), through the development of ad-hoc classes integrated in the standard Geant4 toolkit. This work has been developed in the framework of the PRaVDA (Proton Radiotherapy Verification and Dosimetry Applications) collaboration. The PRaVDA collaboration was established in 2013 to develop a fully solid-state instrument for dosimetry and imaging in the treatment of cancer using proton therapy, and specifically proton Computed Tomography (pCT) [3]. CMOS APSs are integrated in the PRaVDA instrument, along with single-sided silicon strip sensors, in an energyrange detector to measure distal proton residual energies in pCT [4][5][6]. To support the design stage -1 -of the CMOS-based energy-range detector, accurate Monte Carlo simulations, including the charge generation, diffusion and sharing processes, following particle interaction, as well as the full readout chain of the CMOS sensors, were developed.
The remainder of this paper is organised as follows. The CMOS detector used for this work is described in section 2, together with details on the Monte Carlo simulation performed and experiments carried out to validate the charge transport model. Details of the charge transport model proposed are reported in section 3 and, finally, a comparison between experimental and simulated data is reported in section 4 to validate the charge transport model at energies relevant for pCT.

CMOS APS
The wafer-scale APS investigated in this study, named the Dynamic range Adjustable for Medical Imaging Technology (DynAMITe), was manufactured in a 0.18 µm CMOS process by a reticule stitching technique for a total active area of 12.8 cm × 13.1 cm. The DynAMITe pixel array, based on a standard 3T pixel architecture, consists of two different sized diodes meshed in the same pixel matrix. A fine-pitch grid of diodes, offering intrinsic low noise and high spatial resolution, are superimposed to a large-pitch grid of diodes, offering a higher dynamic range. Each cell of the DynAMITe matrix is fitted with multiple diodes: four 0.6 µm diameter photo-diodes placed at the centre of four 50 µm pixels, termed Sub-Pixels, and one 1 µm diameter photo-diode placed at the centre of 100 µm pixels, termed Pixel. The whole matrix comprises 1312 × 1280 Pixels and 2624 × 2560 Sub-Pixels. Both diode matrices can be readout independently or in combination. A schematic representation of the pixel array is shown in figure 1a. For the purposes of this study, only Sub-Pixels are considered. A more detailed description of the pixel architecture, read out modalities and electro-optical performance are reported in [7][8][9].
A cross-sectional view of a Dynamite pixel is shown in figure 1b. For the DynAMITe detector, as in many in standard CMOS technologies, a thin (12 µm) lightly doped p − -type silicon epitaxial layer is grown on a heavily doped p-type substrate (p ++ sub , 735 µm thick). Within the epitaxial layer, n+ wells (n + well ) structures are formed. The p-type epitaxial layer (p − epi ) represents the detector sensitive volume, while the n + well /p − epi diode junction acts a charge collection element. The detector is mainly a field-free volume and only partially depleted across the n + well /p − epi junction ( 1µm), so that charge is collected mainly through a thermal diffusion mechanism. Also, because of the particular doping profile realised across the sensor (p − epi /p ++ sub ), the junction between epitaxial layer and substrate represents a potential barrier, limiting the diffusion of charge generated in the epitaxial layer towards the substrate.

Experiments
Experiments reported here were carried out at the iThemba LABS (SA). The iThemba beam is actively used as a therapeutic facility providing proton beams up to 191 MeV. The maximum beam range achievable at the patient position, or iso-centre, is 240±0.4 mm range with a Full Width at Half Maximum (FWHM) of 25±1.0 mm (measured as 50% of maximum dose on the distal side of the Bragg peak in water). The large area beam (10 cm diameter) is achieved by using a system of -2 - passive scattering components and collimators, while the beam energy can be reduced by graphite attenuators.
The experimental data-set used for this work consists of individual protons measurements, achieved by using a low current (in the range 0.1 to 1 nA, as measured with a transmission ionisation chamber1) and a 10-row Region of Interest (ROI) for the detector readout, resulting in an exposure time of 0.717 ms. The occupancy of the sensor, defined as the ratio of the average number of fired pixels per frame to the number of readout pixels, is in the order of 1%. The pristine beam was reduced in a 30 mm range beam, corresponding to a mean energy of 60 MeV, using a graphite attenuator. This beam was then further reduced by using PMMA attenuators of a thickness of 14, 16, 18, 20 and 22 mm, resulting in a mean proton energy of 38, 35, 30, 26 and 20 MeV, respectively.
Images acquired for these experiments were dark corrected by subtraction of the average of a number of dark frames. Further to this, images were thresholded with respect to a reference value chosen as three times the noise level-as measured soon before the experiment, and equal to 19 DN (corresponding to 1140 e − with a conversion gain of 59 e − /DN). A clustering algorithm was used to account for single hit events spread over multiple pixels.

Simulations
Monte Carlo simulations based on the Geant4 toolkit were conducted. Beam modelling and nozzle geometry were not accounted for in the simulations. The iThemba beam was modelled assuming a 10 cm diameter beam with proton energy (E) extracted from a Gaussian distribution (with a mean 1The exact value of the ratio of proton current to the ionisation chamber beam current is not known, but could be estimated to be in the order of 1%. -3 -energy µ E = 190.8 MeV and standard deviation σ E = 1.53 MeV). Proton emission angle (α) was modelled as extracted from a Gaussian distribution with mean value µ α = 0 and standard deviation σ α = 0.57 deg. The pristine beam was reduced by using a graphite attenuator, whose thickness (124 mm) was adjusted to reproduce the experimental depth-dose curve and range measurement.
Custom classes were added to the standard Geant4 libraries to reproduce the detector response, described in section 3. Signals generated in each detector pixel were recorded from these simulations and the same cluster algorithm, used for the experiment data, was applied to recover the total signal generated in each event.

Charge transport in CMOS APSs
In order to reproduce the detector response to the interaction with radiation, a charge transport model has been developed, consisting of the following steps: 1. Energy deposition in the sensitive volume of the detector is scored by using the standard Geant4 libraries. Individual proton interactions are approximated with a 100-nm step size. Ionisation events, occurring within this step size, result in the generation of a number of secondary electrons, which are tracked as they slow down and their energy deposition recorded.
2. The energy deposited in individual interactions in the sensitive volume is then converted into number of e − /h pairs, assuming an electron-hole pair creation energy for silicon E e−h = 3.6 eV/pair [10].
3. Charge diffusion in the field-free detector occurs, leading to the creation of an electronic cloud.
4. Collection of this electronic cloud is then performed by accounting for collection efficiency as function of the depth of interaction.
5. An artificial process of charge sharing amongst adjacent pixel is introduced, in order to expand the width of the charge cloud.
6. The number of collected electrons are sampled across the detector pixel matrix.
7. The effect of the detector electronics is then included by using the measured detector conversion gain and noise.
In the following sections, each of the steps above will be analysed in detail and a comparison between simulated and measured detector response is shown in section 4, to confirm the validity of this model.

Charge diffusion
From Geant4 standard libraries energy deposition, arising from energy deposition events along a proton track, at a point in space (x 0 , y 0 , z 0 ), where x, y defines the detector plane and z is the orthogonal coordinates to these two, or the beam direction, can be recorded. Energy deposition events, from the primary impinging particle and from all the secondary particles generated, may produce ionisation, generating charge carriers at a mean rate of ≈3.6 eV/electron, equal to three -4 - times the energy band gap of silicon [10]. Ionisation electrons generate a spherical electron cloud with an initial radius σ i given by: with k=0.0062 µm/keV, α = 1.75 and E[keV] being the energy deposited at location (x 0 , y 0 , z 0 ). Equation (3.1) follows directly from the well know energy-range relation R = kE α , where both constants depend on the material in which ionisation is generated [11]. Such generated charge will then undergo thermal diffusion over 4π, as the CMOS APSs are field-free, until it eventually recombines or reaches a collection point. The model used in this work to reproduce charge diffusion in the APS field-free layer was originally developed for field-free regions of CCDs [12][13][14][15]. However, CMOS and CCD sensors present important differences in relation to the ratio between collection electrodes and field-free region. For CCDs, the depleted layer, i.e. the collection electrodes cover the full surface of the pixel, while for most APSs, and -5 -specifically for the device under study, diodes cover only a small fraction of the pixel leading to electrons reaching the top of the pixel in a region not covered by diodes to be reflected back towards the substrate. For this reason, the model used in this work to describe charge diffusion needs to be corrected by a further process described in section 3.3. The charge cloud distribution is approximated in this work by a Gaussian distribution, resulting from charge diffusion in a field-free volume, and its width σ ff is described by the following equation: where z ff is the thickness of the field-free region and z a = z 0 − z ff is the difference between the depth of the interaction z 0 and the thickness of the field-free region. The dependence of σ ff with the depth of interaction z 0 is shown in figure 2a for a 12µm-thick field-free epitaxial layer. The worst case in terms of charge diffusion is represented when energy deposition occurs at the bottom of the epitaxial layer (z 0 = 0), leading to largest possible cloud (σ ff = 6 µm), whereas the best case scenario occurs for charge generated just below the collection area (z 0 = z ff ), where the charge does not suffer charge diffusion (σ ff = 0). The width of the final Gaussian distribution, σ total , is then the sum of the two contributions mentioned above, initial electron cloud generation and charge diffusion, and can be expressed by the following equation:

Charge collection
Using equations (3.1), (3.2) and (3.3), it is possible to build, for each energy deposition event, a quasi-continuous spatial distribution of the ionisation charge undergoing diffusion. However, not all of this charge can be collected by the pixel diodes. In fact, charge diffusing through the detector volume can undergo recombination processes. The probability for charge recombination or collection is dependent on the specific region of the detector volume where charge is generated. Following from figure 1b, three different collection regions can be identified in a CMOS sensor: 1. A heavily doped n + regions near the top of the detector volume, where the high doping concentration and the physical presence of p + wells, for the realisation of the in-pixel transistors, reduces the carrier lifetime and collection is only partial; 2. A lightly doped epitaxial layer, where the low doping concentration leads to a charge carrier lifetime (τ n ), the time available for charge collection before recombination occurs, which is much smaller than the diode collection time (T c ). The collection efficiency in this region approaches 100%.
3. A heavily doped p + +substrate, where the collection efficiency shows an exponential decay due to the short lifetime of charge carriers. Charge collection efficiency decreases for deeper generation points, as recombination is more likely to occur before this charge can reach the collection diode.
The exponential decay in charge collection efficiency for the substrate can be explained by using the extension of Shockley-Ramo theorem as applied to induced charge in semiconductor detectors [16]. The charge transport equation for the excess of minority carrier generated in the p+ substrate can be written as: where D n is the electron diffusion coefficient and τ n is the charge carrier lifetime. Following [17] and using appropriate boundary conditions2 for equation (3.4), the solution of this equation at the collection electrode, as a function of the charge generation point z 0 , can be expressed as [18]: where q(z 0 ) is the charge reaching the collection electrode, N is the charge generated at t = 0 in z 0 , t epi is the thickness of the epitaxial layer, and L n is called electron diffusion length (L n = D n τ n ). The charge collection profile used in this work is shown in figure 2b. A charge collection efficiency of 80% is assumed for the n+ well region ( 1 µm), followed by a full collection efficiency (100%) in the epitaxial layer [17]. Charge collection in the substrate follows the exponential decay of equation (3.5), using a carrier lifetime τ n of 3.5 × 10 −8 s resulting from the doping concentration of the CMOS sensor used in this work. The charge collection profile of figure 2b is consistent with experimental measurements made on similar CMOS technologies [17].

Charge sharing and digitalisation
As discussed in section 3.1, equation (3.2) may lead to underestimate the charge cloud width as the Gaussian approximation for the charge cloud, developed for CCDs, has some limitations when used to model APSs. To correct for this effect an artificial charge sharing model is introduced. The quasi-continuous distributions of charge collected for each ionising event are sampled on an virtual matrix with pixels corresponding to half of the pixel size of the detector (25 µm). In this way each detector pixel (50 µm) is split in four virtual quadrants, each of these containing one of the four in-pixel diodes. This process is schematically shown in figure 2c, where the actual detector matrix is represented by a continuous line, while the virtual sampling grid is shown as a dotted line for the central pixel. Each energy deposition event is then associated to a pixel quadrant, given the calculated spatial distribution of the charge cloud after diffusion and collection (the energy deposition event is represented by a red dot in figure 2c). This charge is then equally split among the four adjacent diodes, around the relevant pixel quadrant (dark grey dot in figure 2c).
Charge collected per ionising event in each pixel quadrant is then re-sampled over the actual detector matrix, i.e. summing up the charge collected for each of the four pixel diodes. Charge collected per pixel has then to be converted into detector Digital Number (DN) and noise has to be added in. In [7] it has been shown how conversion gain and read noise for the DynAMITe detector are log-normally distributed across the pixel array. Signal conversion, from unit of e − to unit of DN, can be performed by randomly extracting per-pixel conversion values from a log-normal distribution with mean µ gain = 59 ± 1e − /DN and standard deviation σ gain = 20.94 ± 0.02e − /DN. Similarly, 2As boundary conditions for equation (3.4), it is assumed that ∆ n = 0 at the edge of the depletion region z = t epi , and at the back of the collection electrode (z = 0). Also, it assumed that the all the charge entering the depletion layer is collected.

Validation
The model for charge transport and detector response simulation, described above, has been validated against experimental measurements acquired with the DynAMITe detector and described in section 2.2.
Detected energy spectra, resulting from signal generated in the detector by individual protons, are shown in figure 3 for experimental measurements and simulations, for the five PMMA attenuators used. Signal distributions, in units of DN, are largely comparable, in terms of width and peak position, for both measurements and simulations. It is to note, however, that measured spectra 26 and 20 MeV proton beams ( figure 3d and e) show the build-up of a small peak at low signal values, which does not appear in the simulations. This could either be due to Random Telegraph Signal (RTS) produced in the detector by radiation damage [19,20] or it could be related to the detection secondary particles, more prominent at low energies.
In order to provide a quantitative verification of the validity of the charge transport model, the distributions of figure 3 have been fitted to a Landau distribution [21], with the most probable value (peak position), corresponding to the mean value of the Bethe-Bloch energy loss. Figure 4 shows the most probable signal, resulting from a fit of the Landau distributions of figure 3, for measurements (red symbols) and simulations (blue curve) as a function of the attenuator thickness, proportional to proton energy. The agreement is reasonable over the whole range investigated, and, even for the data points at 26 and 20 MeV (see figure 3d and e), where some discrepancies are observable in the spectra, measurements and simulations are in agreement within their errors so providing a qualitative verification of the agreement between the charge transport model and experiments.
A further verification of the validity of the proposed charge sharing model is proposed in figure 5. Measured and simulated cluster size, i.e. size of the detected events, is reported in figure 5 for proton in the range 38 to 20 MeV, showing how simulated event size can accurately reproduce measured event size within their errors.

Conclusions
A charge transport model to simulate the response of CMOS APSs to ionising radiation has been developed in the framework of the Geant4 simulation toolkit, in order to support the design stage of a CMOS-based energy-range detector for pCT. The charge transport model included charge generation, diffusion, collection, and digitalisation. The diffusion model used in this work (Gaussian model) was originally developed for CCDs and presented some limitations when used to model the charge cloud width in APSs. To overcome this issue, a non-physically driven artificial charge sharing was introduced. Other authors applied the Gaussian model to APSs and showed the advantages of using a non-physically driven parametric model [22,23]. Other approaches to accurate simulations of charge diffusion in APSs include combining simulation of physics processes (e.g. Geant4) with Technology Computer-Aided Design (TCAD) simulation tools, routinely used for process and device-level design, (e.g. SENTAURUS-TCAD.) [24], or computing the random walk of minority -8 -    . Measured (red symbols) and simulated (blue symbols) cluster size for individual proton events in the DynAMITe detector, when exposed to protons in the range 38-20 MeV.
-10 -carriers [25]. Although very promising, both approaches result time consuming and CPU demanding for many applications. The simulated response of a CMOS APSs to protons in the range 18-38 MeV has been compared with measurements. Simulated and measured data show a good agreement over the whole energy range investigated, both in terms of detected signal and event size, providing a quantitative validation of the charge transport model proposed in this study, which can reliably reproduce the detector response in the energy range of interest for proton CT with sufficient accuracy. Although, the charge transport model presented in this work has been validated with protons in a specific energy range, the model could be seamlessly extended to other radiation fields and energy ranges.