Combining TCAD and Monte Carlo Methods to Simulate CMOS Pixel Sensors with a Small Collection Electrode using the Allpix Squared Framework

Combining electrostatic field simulations with Monte Carlo methods enables realistic modeling of the detector response for novel monolithic silicon detectors with strongly non-linear electric fields. Both the precise field description and the inclusion of Landau fluctuations and production of secondary particles in the sensor are crucial ingredients for the understanding and reproduction of detector characteristics. In this paper, a CMOS pixel sensor with small collection electrode design, implemented in a high-resistivity epitaxial layer, is simulated by integrating a detailed electric field model from finite element TCAD into a Monte Carlo based simulation with the Allpix$^2$ framework. The simulation results are compared to data recorded in test-beam measurements and very good agreement is found for various quantities such as cluster size, spatial resolution and efficiency. Furthermore, the observables are studied as a function of the intra-pixel incidence position to enable a detailed comparison with the detector behavior observed in data. The validation of such simulations is fundamental for modeling the detector response and for predicting the performance of future prototype designs. Moreover, visualization plots extracted from the charge carrier drift model of the framework can aid in understanding the charge propagation behavior in different regions of the sensor.


Introduction
Integrated monolithic CMOS technologies with small collection electrodes [1] are emerging technologies enabling advances in the design of next-generation highperformance silicon vertex and tracking detectors for highenergy physics. These technologies have allowed significant reductions in the material budget with respect to hybrid pixel detectors, while improving the signal-to-noise ratio and the position resolution that is achievable with CMOS sensors.
However, the simulation of such devices remains challenging due to the complex field configuration in the sensor. Advanced simulation tools are required to understand and model the performance of detectors built in these technologies and to optimize the design of future prototypes. This paper presents a simulation performed with a combination of commonly used tools employed in silicon detector simulation. The Allpix 2 framework [2] is used to combine TCAD-simulated electric fields with a Geant4 [3][4][5] simulation of the particle interaction with matter, to investigate the behavior of high-resistivity CMOS detectors and to compare the predicted performance with measurements recorded in a particle beam.
This allows direct access to detector performance parameters such as spatial resolution and detection efficiency by taking into account the stochastic nature of the initial energy deposition. While many of these properties could also be investigated by advanced TCAD transient simulations, this approach is not practical owing to the high computing time for a single event and the high-statistics samples required to evaluate the effects related to the strong variation of the electric field in three dimensions.
Instead, a simplified charge transport algorithm is used, taking as an input the electrostatic field map calculated by the TCAD simulation of the complex field configuration within the sensor. The algorithm takes into account effects like Landau fluctuations in the energy deposition and the production of secondary particles such as delta rays. With event simulation rates of several tens of Hertz, this allows for the generation of high-statistics samples necessary for detailed studies of the detector behavior.
The paper is structured as follows. Section 2 provides a brief overview of the CMOS process under investigation, while the detector properties and the simulated setup are introduced in Section 3. The simulation is described in detail in Section 4, while Section 5 introduces the reconstruction of physical properties from the detector response. The sensitivity of the simulation to a range of parameters is examined in Section 6. The simulation is validated using data recorded in test-beam measurements in Section 7, while performance quantities are derived in Section 8 and compared with the values obtained from data. Finally, Section 9 summarizes the results and provides an outlook for future investigations of this technology.

The High-Resistivity CMOS Process
Monolithic CMOS technologies incorporating the readout electronics in the sensor are attractive candidates for new detector designs to simplify the production and to benefit from a reduction of the material budget. By integrating the CMOS logic in doped wells separated from the inversely doped signal collection electrode, the size of the latter can be minimized as illustrated in Figure 1. The small collection electrode design allows the sensor capacitance to be reduced down to the order of fF, enabling detector designs with low noise and detection thresholds, low analog power consumption, and large signal-to-noise ratio (SNR) [1].
Implemented in a standard silicon substrate, only a small depleted region evolves around the pn-junction surrounding the collection electrode when applying a bias voltage between the doped wells and the backside of the sensor. The applicable bias voltage is limited to −6 V by the process-specific breakdown voltage of the NMOS transistors [7]. In order to achieve a sizable depletion volume around the collection electrode, an epitaxial layer with high resistivity silicon can be used.
The size of the depleted region forming in this epitaxial layer is restricted to the area around the collection electrode and, without additional modifications of the process, no full depletion of the sensor volume is achieved. In the CMOS process under investigation, the depleted region has the shape of a bubble as indicated by the white line in Figure 1, resulting in contributions to the overall detector response from both drift and diffusion of charge carriers. In addition, signal contributions are expected from charge carriers that are created in the highly p-doped backside substrate and subsequently diffuse into the epitaxial layer.

Detector Design under Investigation
The Investigator test-chip is an analog test chip that has been developed within the ALICE ITS upgrade [8]. It has been investigated by the CLICdp collaboration to evaluate this technology in terms of sensor performance focussing on precise measurements of spatial resolution and detection efficiency [6,9]. The digitization of signals is performed off-chip in the data acquisition system using one 65 MHz sampling analog-to-digital converter (ADC) per channel which records the full waveform of all detector channels, once a configurable triggering threshold has been exceeded in any of them [10]. It should be noted that the threshold values for data quoted below represent the offline analysis thresholds applied in addition to the triggering threshold of about 120 e.
The chip has a total thickness of 100 µm. The upper 25 µm of the sensor, below the implants, consist of the epitaxially grown silicon with a resisitiviy of 1 − 8 kΩ cm Sensor¨¨B PCB y in which the depleted region forms, while the additional 75 µm represent the undepleted low-resistivity silicon substrate [7].
While the actual detector contains several sub-matrices with 8 × 8 active pixels each, with different layouts such as altered collection electrode size, only one matrix has been simulated and is compared to data. The pixel cells of the chosen matrix have a pitch of 28 µm × 28 µm and feature the following geometrical parameters: the distance between the p-wells and the collection electrode is 3 µm and an octagonal collection electrode with a size of 2 µm × 2 µm is placed in the center of the pixel cell. A bias voltage of −6 V is applied to the p-wells and a positive voltage of 0.8 V is applied to the collection electrode itself. The simulated detector is placed on a printed circuit board (PCB) as visualized in Figure 2.

Simulation Flow
In the following section, the simulation of the detector in the Allpix 2 framework is described. In order to avoid simulating a full beam telescope setup and performing a track reconstruction, the capabilities of the framework to record the Monte Carlo truth information about primary and secondary particles are exploited.
Consequently, only a single CMOS detector and the source of ionizing particles are simulated as shown in Figure 2. The figure depicts the overlay of many events, as only a single primary particle is simulated in each event.
The following sections describe the individual steps of the simulation in detail, providing information on the configuration of the respective Allpix 2 modules where applicable and relevant.

Electrostatic Field Modeling with TCAD
The electrostatic field in the epitaxial layer of the sensor is modeled using a three-dimensional TCAD simulation. The doping profile is taken from [9,11] and resembles the technology described in Section 2, with the detector geometry introduced in Section 3. The simulation comprises a single pixel cell, and periodic boundary conditions allow the field to be replicated over the entire sensor. Figure 3 shows a visualization of the magnitude of the electric field in the three-dimensional pixel cell, with the corresponding voltages applied to the terminals via metal contacts indicated as gray structures. A low electric field is present in the p-well rings as indicated by the blue region on the surface of the simulated pixel cell. The center of the p-well rings is fitted with a squared opening that contains the collection electrode with a high-field region evolving around it.
The strong inhomogeneities of the electric field in different regions of the pixel cell are best observed in a cut through the collection electrode, perpendicular to the sensor surface, as depicted in Figure 4. The high electric field strength close to the pn-junction around the collection electrode decreases rapidly towards the sensor backside and the pixel corners. The white line indicates the depleted volume of the pixel cell. The electric field lines, indicated as black arrows, provide a first insight into the complexity of the field configuration in the sensor and the drift effects induced by this strong non-linearity. The low electric field regions in the pixel corners result in a slower charge carrier drift and an increased impact of diffusion, leading to an enhanced charge sharing which improves the position resolution without the need to reduce the pixel pitch. In the low-resistivity substrate, recombination of charge carriers is a relevant process owing to the higher doping concentration.
The electrostatic field obtained from TCAD is converted to a regularly spaced mesh using the Mesh Converter tool provided with the Allpix 2 framework. This conversion speeds up the look-up of field values during the simulation by several orders of magnitude since all necessary interpolations are already performed offline prior to the simulation. A regular mesh granularity of 0.1 µm is chosen to correctly replicate the field in the high-density mesh regions of the TCAD simulation close to the implant.
It has been verified that the selected granularity correctly replicates the TCAD simulation by comparing the two fields. Using an even finer granularity has not shown any significant improvement on the simulation results. Loading highly granular electrostatic fields in Allpix 2 does not impact the performance of the simulation, but only the memory footprint of the program during execution.

Energy Deposition with Geant4
Allpix 2 provides the DepositionGeant4 module, an interface to Geant4 [3][4][5] which facilitates the simulation of energy deposition in the sensor. A 120 GeV beam of π + incident on the pixel detector is simulated, replicating the beam conditions of the test-beam measurements. The beam points along the positive z -axis, perpendicular to the xy-plane of the detector. The cross section of the beam is chosen to be significantly smaller than the detector surface to suppress effects stemming from the altered charge sharing behavior at the sensor edge. The energy deposited in the sensor by Geant4 is translated into charge carriers with a conversion factor of 3.64 eV per electron-hole pair [12].
The framework also stores the Monte Carlo truth information including secondary particles such as delta rays and their relation to the primary particles. This information can be exploited to establish a link between the incident particles and the electron-hole pairs created in the detector.
The simulation is performed with the Photo-Absorption Ionization model (PAI) [13] to improve the description of energy deposition in thin sensors. This is of importance in spite of the total sensor thickness of 100 µm, since a majority of the charge carriers forming the final signal will originate from the 25 µm epitaxial layer.
Listing 1: Configuration section for the DepositionGeant4 module setting up the particle source in Geant4 for the initial energy deposition in the sensor.
The module configuration used for the Allpix 2 framework is provided in Listing 1.

Charge Carrier Transport
The signal formation is simulated using a simplified model for charge carrier transport based on the notion of collected charges. The electron-hole pairs created by the incident particle are propagated along the electric field lines through the sensor volume using an adaptive fourthorder Runge-Kutta-Fehlberg (RKF) method [14] and a mobility parametrization which depends on the electric field vector [15]. The RKF method adapts the simulated time step depending on the position uncertainty derived from a fifth-order error estimation; the allowed range for time steps was set to 0.5 ps ≤ ∆t ≤ 0.5 ns.
While this model is not expected to reproduce a realistic time dependence of the signal, the final state of charge collected at the sensor implants is equivalent to the integrated induced current over the respective drift time. This approximation is valid since the Shockley-Ramo weighting field [16,17] is negligible in most of the sensor volume owing to the small ratio between signal collection electrode size and sensor thickness.
In the upper 25 µm of the sensor the charge carrier motion is a superposition of drift and diffusion, while in the lower 75 µm the charge carriers are only subject to random diffusion as the electric field is negligible.
The propagation algorithm is halted after 22.5 ns, the so-called integration time, and all charge carriers within a volume of 3 × 3 × 2 µm 3 around each of the signal collection electrodes are attributed to the respective pixel signal. The volume has been chosen to cover the electrode implant itself as well as an additional volume accounting for the uncertainty in the final position of the transported charge carriers. The integration time has been chosen such that the simulation produces clusters with the same most probable value (MPV) for the cluster charge as obtained from data. This aims to emulate the physical process of charge carrier recombination in the silicon substrate, which might be modeled directly in future simulations as briefly discussed in Section 9. The systematic uncertainty introduced by this approach is discussed in Section 6.
Charge carriers are transported in groups of five instead of individually to speed up the simulation process. The group size has been chosen such that an adequate number of transport steps is retained with the expected MPV for the signal of around 1.5 ke. It has been verified that this simplification does not affect the simulation result as further elaborated in Section 6. Figure 5 visualizes this transport model and shows the collection of charge carriers at the electrodes of the sensor. In this representation, only electrons that have reached a sensor implant within the integration time are shown. Electrons that are still in motion as well as holes are suppressed. The motion of each group of charge carriers is represented by one line and is shown at different integration times after the initial energy deposition. Here, the incident particle traversed the detector along the z-axis through the center of one pixel cell.
After the first few hundred picoseconds, only charge carriers in the vicinity of the electrode are collected. The straight lines indicate that their motion is dominated by drift. With increasing integration time, the motion patterns of further groups of charge carriers arriving at the implant exhibit a strong contribution from diffusion as indicated by the numerous kinks in the respective paths. After about 15 ns, lateral motion enables some charge carriers to be collected in the two adjacent pixel cells.
The line graphs also allow visual distinction between the substrate and the epitaxially grown high-resistivity layer, which ends about 25 µm from the top of the sensor. A faster drift motion can be observed in the high-field region close to the backside of the epitaxial layer as straight lines; the contribution from substrate charge carriers diffusing into the epitaxial layer starts only after approximately 10 ns.
In Figure 6, a three-dimensional representation of the line plot at 20 ns is presented. The lines end at five different points, each representing a different collection electrode.
The configuration provided in Listing 2 has been used for the charge carrier transport. Settings for creating line graphs of the charge carrier motion can be found in the Allpix 2 user manual available from the project website [18].

Digitization of Signals
To simulate the response of the readout electronics, the charge carriers accumulated in the region around the signal collection electrode during the integration time are transformed into a digital signal. While the detector under investigation uses off-chip ADCs for the signal digitization as described in Section 3, the simulation aims to simulate an on-chip per-pixel threshold using the DefaultDigitizer module of Allpix 2 . Equivalent noise values have been used where applicable, as discussed below.
[DefaultDigitizer] electronics_noise = 10e threshold = 40e threshold_smearing = 5e An additional signal contribution, randomly drawn from a Gaussian distribution with a width of 10 e and a mean of 0 e is added to the signal to account for electronics noise present during digitization. The applied threshold is varied between 40 e and 700 e, and a threshold dispersion, sampled from a Gaussian distribution with a width of 5 e and a mean of 0 e, is added. For simplicity, the threshold dispersion is not a fixed offset calculated per-pixel, but randomly chosen per pixel hit. The setup of the module is summarized in Listing 3.

Data Processing and Storage
The simulation results are stored in ROOT [19] trees using the ROOTObjectWriter module. In order to speed up the process, the simulation is performed in two separate steps. In the first step, the energy deposition, charge carrier transport and summing of charge at the collection electrodes is performed. The result of this step is stored to disk.  In a second step, the ROOTObjectReader is used to read the information from the respective file and the final digitization step is performed. This allows to re-run this final section of the simulation on the full set of Monte Carlo events with different settings applied without the need to recompute the drift motions. A full threshold scan, performed on the same set of initial simulation events, thus only takes a couple of minutes instead of several hours required to create the initial data set. Since the threshold scan performed on the test-beam data has also been performed offline on the same data set [9], this is an adequate simplification of the simulation.
The central simulation data set comprises about 2.5 million primary events which have been reprocessed for every threshold setting. In addition, several smaller data sets with different integration times have been produced in order to optimize agreement with data as discussed in Section 6.

Reconstruction and Analysis
In the following, the reconstruction and analysis of the Monte Carlo events are discussed. The simulation was set up using known, independent parameters of the measurement setup, such as track resolution or charge threshold. Only the cluster charge MPV was used as direct observable provided by the detector to tune the simulation. All parameters were fixed before comparison with data for observables used to quantify the performance, such as cluster size, position resolution and efficiency. This blinded approach avoids drawing premature conclusions from the figures of merit and thus influencing the parameter optimization. Using only the MPV of the cluster charge for calibrating the simulation minimizes the correlation between simulation and data, and maximizes the prediction power of the simulation.

Reference tracks
The Monte Carlo truth information provided by the Allpix 2 framework is used as reference track information.
All registered particles in the sensor are filtered and only primary particles entering the sensor from the outside, i.e. those without a parent particle, are selected for further analysis. This set of particles represents external tracks, and their position in the mid-plane of the sensor is calculated by linearly interpolating their entry and exit points registered by the framework. This position is then convolved with the track resolution at the device under test (DUT) of 2.0 µm, in accordance with the value obtained for the beam telescope used for the acquisition of the testbeam data [20].

Clustering
The pixel hits registered by the simulated detector are grouped into clusters by starting with a seed pixel and adding all directly adjacent pixel hits to the cluster until no additional neighbors are found. This already allows basic properties of the simulation to be compared with data, namely cluster size as well as the shape of the cluster charge distribution.
The total cluster charge is given by the sum of the individual pixel charges of the cluster. Its comparison with data allows the required integration time in the simplified simulation model to be adjusted to achieve the same integrated charge as seen in data. This procedure is described in detail in Section 6.
The cluster size is defined as the total number of pixels contained in the respective cluster. It has a strong dependence on the drift and diffusion of the charge carriers in the sensors and is the primary measure for charge sharing between pixel cells. It thus allows evaluation of the performance of the simulation, e.g. how well the electric field is modeled.

Reconstruction of the Cluster Position
For assessing the performance of the detector, a particle incidence position has to be extracted from the cluster information available. To replicate the analysis performed for the test-beam data, the charge-weighted centerof-gravity position of the cluster is corrected for non-linear charge sharing by an η algorithm [21].
Since the η distribution represents the charge sharing between two pixels only, for each cluster the two pixels with the highest charge, Q 1 and Q 2 , are chosen to construct the η variable independently in x and y: where k i is the relative position between the two pixel centers. An example of the η distribution in x is depicted in Figure 7 for a pixel charge threshold of 40 e.

Systematic Uncertainties
The sensitivity of the simulation to different input parameters has been examined by varying the values within their respective uncertainties, if known, or within a reasonable range otherwise. The impact on the reconstructed observables was investigated. While some parameters exhibit little or no effect on the simulation results, others have a strong influence on the outcome.

Free parameters
For the initial deposition of energy in the sensor, the influence of the maximum allowed step length of tracking primary and secondary particles through the sensor material has been evaluated by varying the respective value between 0.1 µm and 5 µm, and no significant difference was observed. Since large parts of the sensor volume are undepleted, a strong impact of diffusion is expected which smears the initial position of the charge carriers.
The charge carrier transport is mainly dominated by the precision of the numeric integration and its granularity. The number of charge carriers transported as group has been varied from a single charge carrier up to ten per group in order to study possible effects on the distribution at the implants. The effect on the reconstruction observables is found to be negligible.

Parameters constrained by measurements
The behavior of the sensor has been shown to be very sensitive to the simulated physical properties of the CMOS sensor, i.e. the thickness of the epitaxial layer as well as the modeled electric field. Even small changes in the sensor design, such as a more simplistic approximation of the implant doping profiles in the TCAD design cause large changes in the resulting cluster size distributions and position resolution. It is therefore of paramount importance to model the sensor as precisely as possible and to constrain the different parameters in TCAD by additional measurements [22]. The low-field regions found in the corners of the pixel cell visible in Figures 3 and 4 are strongly influenced by these modifications, and their contribution to the detector signal changes accordingly. The integration time currently used to stop the transport of charge carriers is also linked to the sensor design, since it is used to emulate an effective recombination of charge carriers. Their lifetime in the different regions of the sensor is dominated by the respective doping concentration, and potentially affected by the silicon wafer production process. Since this was not modeled in detail for this simulation, the integration time was chosen such that the MPV of the cluster charge matched the value obtained from data as discussed in Section 4. The corresponding uncertainty on the charge calibration of the reference data has therefore to be taken into account as systematic uncertainty of the simulation by comparing the cluster charge MPV for different integration times to the value obtained from data as shown in Figure 8. Here, the hatched band represents an assumed uncertainty of ±50 e on the charge calibration of data [7,9]. This translates to an uncertainty on the integration time of 22.5 +1.5 −1.3 ns, which is propagated as systematic uncertainty to the results presented in this paper.
It has been observed that the overall agreement between data and simulation seems to improve for lower integration times, which might indicate either an offset in the absolute charge calibration of data or an insufficient modeling of the signal formation processes in silicon.
Also the charge threshold applied to the individual pixels has a strong impact on both the cluster size and the intrinsic resolution, with decreasing influence towards higher thresholds. At a threshold of 40 e, a change of as little as ±5 e visibly alters the observables. Since the absolute charge calibration and the threshold value in electrons are fully correlated, the uncertainty on the applied threshold has been taken into account by varying the two parameters simultaneously and by calculating the total uncertainty arising from the variations.
A variation of the threshold dispersion and electronics noise of up to 10 e at a threshold of 40 e yielded no observable effect. The values for noise and threshold dispersion have been estimated from the evaluation of the full waveform in data [9].
The residual width and the final intrinsic resolution depend on the resolution of the reference tracks at the position of the detector under investigation. This resolution has been determined for the test-beam data used, and a variation of ±0.2 µm around this value shifts the obtained resolution accordingly. This strong influence arises from the fact that the two values are of similar size.
In summary, while the free parameters of the simulation have little to no influence on the final result when varied within a reasonable range, several parameters show a high sensitivity but are constrained by measurements.

Validation With Test-Beam Data
The simulation is compared to data recorded with the Investigator chip, described in Section 3, at the CERN SPS accelerator with a 120 GeV π + beam. A total of 25660 tracks through the region of interest have been recorded, mainly limited by the very small active area of the DUT and the dead time of the data acquisition system used. More details about the test-beam setup, data samples and the analysis of data used for comparison in this paper can be found in [6,9].

Cluster Charge
The cluster charge distributions for both simulation and data at a charge threshold of 120 e are shown in Figure 9. The distributions are fitted with the convolution of a Gaussian and Landau function. The MPV is 1.42 ke for both data and simulation, and the width of the Gaussian is 0.21 ke/0.22 ke for data/simulation, respectively. A good agreement between data and simulation is observed, as also indicated by the ratio of the two distributions displayed in the lower part of the figure. While the MPV has been tuned to match data using the integration time of the simulation as discussed in Section 4, the agreement of the shapes indicates that the energy deposition in the relevant parts of the sensor as well as the collection of charge carriers is well-modeled by the simulation. The data distribution exhibits some fluctuations owing to the low statistics of the sample.

Cluster Size
The distribution of the total cluster size at a threshold of 120 e for simulation and experiment is presented in Figure 10. Qualitatively, the distributions are in good agreement. A possible source of the observed deviations for individual cluster sizes are uncertainties in the modeled electric field of the sensor as discussed in Section 6.
The projection of the cluster size in x and y, depicted in Figure 11, provides additional details about the charge sharing process. Data and simulation agree well, but a small difference between the distributions in x and y can be observed in data despite the symmetry of the pixel cell layout. It has been verified that this does not stem from a remaining misalignment in data by repeating the simulation with a sensor rotated around the x axis by up to ±15°i n an attempt to reproduce the difference. The deviation might be a result of the non-symmetric layout of the circuitry in the Investigator pixel. While the p-well structure has been designed to be fully symmetric in x and y, the layout of the circuitry placed in the p-wells is not symmetric, which is a possible source of the asymmetry.
The cluster size distribution is a precise measure for charge sharing as confirmed by the intra-pixel representation of the total cluster size presented in Figure 12. For the simulation, the Monte Carlo truth information is exploited to produce a multi-pixel map indicating the mean cluster size as a function of the particle incidence position within the pixel cells. Likewise, the reference track supplied by the beam telescope is used to obtain the particle incidence position for data. To increase statistics, data events from the full active matrix are folded into a single pixel cell, which is displayed in the upper-right quarter of Figure 12.
The largest clusters originate from the pixel corners since the low electric field between pixel implants results in a strong contribution from diffusion of charge carriers. Single-pixel clusters, on the other hand, are almost exclusively produced if the incident particle traverses the sensor very close to the center of the pixel cell.
While the overall mean size distribution is faithfully reproduced in the simulation, minor discrepancies in the pixel corners are visible. The transition from four to threepixel clusters represented by the yellow regions is more apparent in simulation than in data. The same holds true for the transition between two to three pixel clusters corresponding to the turquoise regions in Figure 12. Particles penetrating the sensor at the corners of a pixel cell, for example, are more likely to give rise to clusters with size four in data compared to simulation. This observation is in line with the higher number of clusters with size four in the cluster size distribution displayed in Figure 10. Moreover, the cluster size is particularly sensitive to a mis-modeling in the pixel corners as the diffusion of charge carriers to neighboring pixel cells is most likely if the incident particle enters the sensor at the corners between four cells. Most notably, small modifications in the electric field in the pixel corners are capable of inhibiting or enhancing the motion of charge carriers to neighboring cells causing deviations in cluster size by up to two units as there are two cells directly adjacent to one corner. As discussed in the previous section, the low field regions in the pixel corners are strongly influenced by the exact doping profile of the sensor.  The mean cluster size has been studied as a function of the applied charge threshold. Figure 13 shows the curves for data and simulation. In addition, a simulation with a linear electric field replacing the TCAD model in the epitaxial layer is plotted as a dashed line for comparison. By increasing the threshold, the mean cluster size shifts to smaller values as individual pixels fall below the charge threshold. Data and simulation match well down to very low thresholds, with a maximum deviation of about 6 % at very low thresholds, while the simulation with a linear electric field produces incompatible results. This deviation from the experimental results demonstrates the significance of a precise modeling of the electric field for this type of detector. Similar results have been obtained for the mean projected cluster sizes along the x and y coordinates. Figure 14 displays 2 × 2 pixel maps of the mean projected cluster size in x and y as a function of the particle incidence position at a threshold of 40 e. Instead of the uniform bands along the respective coordinate expected for uncorrelated observables, eye-shaped structures reveal a correlation between charge sharing along the two dimensions caused by the inhomogeneous electric field and the bubble-shaped depletion region described in Section 2. The same effect is observed in data as demonstrated in [9]. With increasing threshold, charge sharing effects are suppressed and the correlation between the mean cluster size in x and y vanishes.

Detector Performance
Using the reconstructed cluster position and the Monte Carlo truth information from the primary particle, the performance of the CMOS detector is assessed in terms of spatial resolution and hit detection efficiency. The results obtained from simulation are compared to data. Figure 15 shows the residual in x, defined as the difference between the impact point of the incident particle obtained from the Monte Carlo truth and the reconstructed cluster position. The width of the residual is obtained as the root mean square (RMS) of the distribution, evaluated for the central 99.73 % of the histogram, equivalent to ±3σ of a Gaussian distribution, to match the definition used in the data analysis. This allows the width of the distribution to be quantified independently from its shape while providing a statistically robust metric.

Intrinsic Resolution
The spatial resolution is then calculated by quadratically subtracting the track resolution from the residual width, i.e.
The statistical uncertainty on the resolution is calculated using pseudo-experiments. The number of entries in each bin of the residual distribution under consideration is smeared with a Poisson distribution with a mean equivalent to the original bin content. The width obtained from the smeared histogram is stored, and the pseudoexperiment repeated 10 000 times. The statistical uncertainty on the residual width is then taken as the width of the resulting distribution and is propagated to the intrinsic resolution.
The resolution has been studied as a function of the charge threshold applied, shown in Figure 16 for the x and y coordinates separately. With increasing threshold, the information from pixels not passing the threshold is lost, leading to a deterioration of the position resolution. The comparison of data with simulation shows a very good agreement down to a threshold of about 150 e. The discrepancy at lower thresholds is most likely to be a consequence of non-Gaussian noise in the data recorded with the analog prototype chip as well as a result of the simplification of charge carrier lifetimes described in Section 4. The disagreement is of limited importance for practical purposes since a fully integrated sensor is likely to be operated at thresholds above 150 e.
The dashed gray line in Figure 16 again represents a simulation using a linear electric field as approximation, and the deviation from data suggests that this simplification leads to an inadequate description of the CMOS sensor response.

Efficiency
The efficiency of the detector is defined as the number of incident primary particles that can be matched to a reconstructed cluster divided by the total number of primary particles penetrating the detector. A match between an incident particle and a reconstructed cluster is made, if the cluster is located within a radius of 100 µm around the impact point of the incident particle, using the same matching criterion as applied to data.
The statistical uncertainty of the efficiency has been calculated by considering a counting experiment with two possible outcomes: either a matched or an unmatched primary particle track. This results in an uncertainty of where p is the probability of a matched track while N is the total number of experiments conducted. The efficiency obtained from simulation as a function of the particle impact position within a single pixel cell is displayed in Figure 17 for three different thresholds. −0.23 (syst) %. The statistical uncertainty is of the order of 1 × 10 −8 . The remaining inefficiencies are evenly distributed throughout the pixel cell and arise from delta rays which pull the cluster center far away from the particle incidence point. With increasing threshold, inefficiencies start to develop in the pixel corners, as these are the regions with the strongest charge sharing and the largest mean cluster size. The overall hit detection efficiency at the threshold of 450 e shown in Figure 17 (center) decreases to about 97.62 +0. 13 −0.58 (syst) %. At the threshold of 700 e, depicted in Figure 17 (right), a pronounced inefficiency is observed, extending from the pixel corners into the pixel cell and leading to an overall efficiency of 85.96 +0.53 −1.02 (syst) %.
This decrease of efficiency can best be observed as a function of the charge threshold applied, as shown in Figure 18. While the shape of the curve observed in data is reproduced well, a constant offset to the measured values can be observed. This difference can be attributed to fluctuations of the pedestal as well as inefficiencies in the data acquisition system which are not modeled in simulation. The simulation using the linear electric field approximation is found to not correctly model the behavior observed at high threshold values.

Summary & Outlook
In this paper, a combined 3D TCAD and Monte Carlo simulation of a CMOS pixel sensor with small collection electrode design, implemented in a high-resistivity epitaxial layer, has been presented. The simulation combines the results of a three-dimensional electrostatic TCAD simulation with the stochastic description of energy deposition by Geant4 using the Allpix 2 framework. Visualizations of the charge carrier motion in the sensor produced by the simulation framework have been found to be helpful to qualitatively understand the sensor response.
The simulation results have been compared to measurements of a reference detector, recorded in a test-beam, and very good agreement has been observed after tuning the simulation to match the most probable value of the cluster charge measured in data. The simplified charge transport model implemented in Allpix 2 has been shown to be sufficiently precise to replicate the detector performance figures of merit such as efficiency and intrinsic resolution measured in data.
The implemented simulation setup for CMOS sensors will be used for further studies of similar detector prototypes and designs, including different sensor geometries and modified production processes aiming at a full lateral depletion of the epitaxial layer.
In future versions of the Allpix 2 framework, a simulation of charge carrier recombination might be implemented, calculating the lifetime from the respective doping concentration as a function of their position within the sensor. This would allow for an even more realistic description of the charge transport process and would remove the necessity of setting and tuning the integration time for underdepleted detectors.
Furthermore, the simulation could be extended to the detector performance in the timing domain by simulating the charge transport taking into account induced currents using the Shockley-Ramo theorem as possible with the latest version of the Allpix 2 framework.
The presented combination of precise electric field modeling in TCAD and inclusion of statistical fluctuations is also interesting for the simulation of other silicon detector technologies with complex field configurations such as 3D sensors or LGADs.