Simulation of semiconductor detectors in 3D with SolidStateDetectors.jl

The open-source software package SolidStateDetectors.jl to calculate the fields and simulate the drifts of charge carriers in solid state detectors, together with the corresponding pulses, is introduced. The package can perform all calculations in full 3D while it can also make use of detector symmetries. The effect of the surroundings of a detector can also be studied. The package is programmed in the user friendly and performance oriented language Julia, such that 3D field calculations and drift simulations can be executed efficiently and in parallel. While all kinds of semiconductor devices can be simulated, special emphasis is put on germanium detectors. The verification of the package is shown for an n-type segmented point-contact germanium detector. Additional features of SolidStateDetectors.jl planned for the near future are listed.


Introduction
Semiconductor detectors are used widely in industry and in basic as well as applied science. They play an important role in various nuclear and particle physics experiments. A special case are highpurity germanium (HPGe) detectors, for which SolidStateDetectors.jl was originally developed. Due to the extremely low impurity densities achieved for germanium, typically around 10 10 per cubic centimeter, these detectors can have volumes of several hundred cubic centimeters. Germanium detectors are known for their good energy resolution, which is in many cases exploited by only reading out one of two electrodes, often called "core". Cylindrical geometries are common where the core electrode is established either in a central bore hole or as a circular patch on one of the end-plates. In order to minimize the capacitance of such a patch, it is often implemented as small as possible and is called point-contact. Depending on the manufacturer, such detectors are also called "Broad Energy Germanium" (BEGe) detectors.
In AGATA [1] and in GRETINA/GRETA [2,3], arrays of segmented cylindrical HPGe detectors are used for spectroscopic measurements using gamma-ray tracking. The outer surface, also called mantle, of an AGATA detector, which is 10 cm long and 8 cm in diameter, is divided into 36 electrodes and for a germanium detector, this is considered highly segmented. Low-background experiments searching for neutrinoless double beta (0 ) decay [4,5] or dark matter [6][7][8] as well as experiments measuring coherent elastic neutrino-nucleus scattering [9][10][11] exploit the radio-purity, the excellent energy resolution and the low energy threshold of point-contact HPGe detectors. These detectors are usually not segmented and the point-contact is the only electrode which is read out.
In modern experiments, not only the energies deposited in the detectors are recorded but also the time development of the charges induced in the read out electrodes, the signals which form pulses. Pulse shape analysis plays a central role to facilitate position reconstruction in segmented detectors and background rejection in point-contact detectors. In particular, the rejection of events with more than one energy deposition, i.e. multi-site events, in searches for 0 decay reduces the background induced by gamma-rays by about an order of magnitude [12], while keeping a high efficiency for single-site events like 0 decays. In the case of point-contact or BEGe detectors, both G [4,12] and M [5,13] apply pulse shape discrimination based on the so called / parameter, where is the maximum of the current amplitude and is the deposited energy. In an alternative approach, artificial neural networks have been implemented to reject multi-site events recorded in coaxial [12,14] and BEGe [15] detectors. Based on the success of G [16] and M [17] in achieving good sensitivity to 0 decay, the next generation experiment LEGEND [18] will, in its first phase, operate about 200 kg of HPGe detectors under low-background conditions. In its second phase, LEGEND will be upgraded to 1000 kg.
A precise simulation of the detector pulses is essential to fully understand the effects of pulseshape discrimination techniques and their efficiencies, and to provide clean training samples for neural networks. One of the technical applications of simulation is that, during the production of detectors for LEGEND, the cutting of the available crystals can be optimized to make the production of detectors with maximal masses and still reasonably small bias voltages possible. This optimization has to take the active impurity profiles of the crystals into account. Thus, a simulation package is needed, where not only the geometry but also any kind of impurity profile can easily be implemented. In addition, the electric field calculation including the handling of undepleted regions has to be fast, such that many scenarios can be studied in a reasonable amount of time. However, the package is not an engineering tool. It does not provide any checks whether a detector is technically feasible. The package provides the functionality to predict the response of a given detector configuration in a given environment.
Previously published packages to simulate the pulses of HPGe detectors do not take the detector environment into account. In addition, they are either not open source [19,20] or difficult to extend [21,22]. The new open-source pulse-shape simulation package, SolidStateDetectors.jl, fulfills all requirements of LEGEND related research activities and provides a modular software environment that can easily be extended by the user for other applications like different HPGe detector geometries or special silicon detectors.

Detector geometries and simulation procedure
The user defines the "world", including the geometry of the detector and, optionally, its surroundings as well as the electric boundary conditions in structured text configuration files. Constructive solid geometry (CSG) is used to define different shapes and to construct all objects. The coordinate system, either cylindrical or Cartesian, is also defined in the configuration file as well as the boundaries of the world. SolidStateDetectors.jl provides example configuration files for selected common HPGe detectors, including a segmented true coaxial [23], a segmented BEGe [24] and an inverted coaxial [25] detector . It also provides an interface to read in configuration files from M Siggen [22]. The simulation can be divided into two main parts: 1) the calculation of the static electric properties of a given experimental setup and 2) the drift of charge carriers inducing signals on the electrodes. Both of these are implemented in the programming language ulıa in a modular manner. The ulıa language has a high-level syntax leading to an easy to use environment but it is still designed for high performance [26].
The electric field is calculated once for a given detector geometry at the beginning of the simulation, based on the electric properties of the system. The main properties are the impurity density profile of the crystal and the fixed potentials of the electrodes. In addition, the user can define volumes with fixed charges on the detectors and surrounding materials at fixed electrical potentials, like grounded holding structures. Similarly, the weighting potentials, which are used to determine the signals induced on the electrodes by the charge carriers, are calculated once for each detector configuration.
For an event-by-event simulation, the interaction of radiation with the detector and surrounding material can be simulated with a dedicated software package like G 4 [27]. The positions of the individual energy depositions (hits) of each event are used as input to the pulse generation. SolidStateDetectors.jl provides flexible clustering of the input hits through the service package Clustering.jl; this can be performed when the hits are read in. For each hit or cluster, the induced charge carriers are drifted in user defined time steps, default is 1 ns, according to the previously calculated electric field and the implemented drift velocity model. At each step, the induced signals on the electrodes are computed using the weighting potentials.

Electric field and weighting potentials
The electric field, E(r), is calculated numerically as the first derivative of the electric potential, Φ(r), which itself is calculated numerically according to Gauss' law: where 0 is the vacuum permittivity, (r) and (r) describe the relative permittivity and the charge density as functions of the position r. The input on (r), (r) and the boundary conditions is specified by the user in the configuration file. The charge density inside a detector has two different components: (r) = imp (r) + fix (r), where imp (r) is due to the impurity density of the semiconductor and fix (r) is an optional fixed contribution, e.g. a charged surface layer. The numerical calculation is performed using a successive over-relaxation (SOR) algorithm on an adaptive grid with red-black division for parallelization. The two nearest neighbors in each dimension are used to numerically update the potential, Φ , on the grid point . At the boundary of the world, a prescription how to update the potential of the point , which is located on the boundary, The latter is the baseline design for the LEGEND experiment because such detectors have excellent pulse-shape discrimination properties and large masses. is needed. An extra point +1 is introduced, for which the user can choose for each axis individually the prescription: reflecting (Φ +1 = Φ −1 ), periodic (Φ +1 = Φ 1 ), fixed (Φ +1 = constant, normally Φ +1 = 0), decaying (Φ +1 = Φ / +1 , where ( +1 ) is the distance between point and the origin in this dimension). Reflective or periodic are good choices if the symmetries of a system are to be used. Periodic is especially useful for the -direction, where it allows to calculate the fields of rotationally symmetric detectors in 2D instead of in full 3D.
At the beginning of the calculation, a very coarse grid is initiated. Once the SOR algorithm has converged, i.e. changes are below a user defined value, grid points are inserted in a dimension if the difference in potential for two neighboring points is above a user defined threshold . The default number of such grid refinements is three. However, the user can set the number to a higher value to ensure that the targeted precision given by the threshold is achieved. The gradual adaptive refinement of the grid avoids large gradients between grid points and ensures fast convergence of the calculation as it does not create a large number of points in volumes of moderately changing potentials. The electric field at a given position between grid points is evaluated by linear interpolation.
There is an option that during an iteration within the SOR algorithm, grid points are marked as undepleted. As the net charge carrier density is zero in the undepleted regions of a semiconductor, imp (r) is set to zero for such grid points. This feature allows studying the development of the depletion zone with increasing bias voltage for different detector geometries and impurity distributions in order to optimize detector configurations. Especially point-contact detectors can, depending on the impurity density profile, develop volumes which cannot be depleted. After a crystal is produced, it is essential to cut it such that this problem is avoided. Simulation is the best guide for this.
The weighting potential, Φ (r), determines the size of the signal induced on the electrode as a function of r. It is calculated by solving with the boundary conditions Φ (r ) = 1 and Φ (r ) = 0 for all r on electrode and all r on electrodes ≠ . The same SOR algorithm as used to calculate the electric potential is used to calculate the weighting potentials. Figure 1 depicts the n-type segmented BEGe detector [24] which is used as an example throughout this paper. This detector features a point contact (core) at the top end-plate. The core electrode is surrounded by a ring covered with a passivation layer. It has a rather unusual segmentation, which was chosen to study the signal development in BEGe detectors and makes it particularly useful to verify a simulation package. This segmentation is four-fold: Three equal segments (1=red, 2=yellow, 3=grey) cover about one sixth of the surface each. The remaining surface is covered by one large segment (4=magenta), which is closed at the bottom. The center of the bottom plate is the origin of a cylindrical coordinate system with the -axis pointing towards the core contact. The left edge of segment 1, looking from the side, defines = 0 • .
The electric potential and the electric field of the detector as simulated with SolidStateDetectors.jl are shown in Figs. 2a) and b) for the situation that the detector is surrounded by vacuum. The high electric field underneath the core contact is typical for BEGe detectors.  An infinitely long true-coaxial geometry was implemented to validate the numerical field calculation by comparing the result obtained with SolidStateDetectors.jl to the analytical solution available for this configuration. The inner (outer) radius of the germanium detector was set to 5 mm (35 mm) and the charge density (r) was defined as proportional to 2 inside the crystal and zero elsewhere. The potential on the inner (outer) mantle was set to 0 V (10 V). The numerical calculation was performed using both cylindrical and Cartesian coordinates with boundary conditions reflecting in and periodic in . The symmetries of the detector were used. After the default of three refinements, the mean distance between grid points became ( : A comparison of numerically and analytically calculated values of the electric potential is presented in Fig. 3. The root mean square, , difference between the numerical and analytical solutions is very small for the calculation performed in Cartesian coordinates, = 0.025 V, and almost zero for cylindrical coordinates, = 2.1 · 10 −6 V. These results were obtained with the automated termination of iterations and grid refinements. In both cases, the could be further reduced if iterations and grid refinements beyond the defaults were initiated. However, the precision obtained with the defaults would be sufficient for any real-life detector. The capacitance, , of a detector can be calculated as = 2 / 2 , where is the applied bias voltage and is the energy stored in the generated electric field: The integral becomes a sum over all grid points of the world. The procedure is verified by comparing numerically to analytically calculated capacitance values for the infinitely long true-coaxial detector and, in addition, for an infinitely long and wide parallel plate detector. The ratios of the numerical (cyl = cylindrical coordinates, car = Cartesian coordinates) and analytical (true) values are coax cyl / coax true = 0.9997, coax car / coax true = 1.01, plate car / plate true = 1.006. These ratios asymptotically approach unity if the number of iterations and grid refinements is increased beyond the default values. The comparisons described here are part of the automatic tests, which are performed every time the source code of SolidStateDetectors.jl is changed.
In SolidStateDetectors.jl, the relative permittivity (r) is a position dependent parameter. This allows the user to implement the surroundings of the detector, e.g. infrared shields or support structures in close proximity to the crystal. While HPGe detectors are typically operated in vacuum cryostats, collaborations like G [4] or LEGEND submerge their detectors directly in liquid argon. Figure 4 demonstrates the importance of taking the surroundings of a detector into account. Especially, volumes close to the passivated surfaces exhibit non-negligible differences in the electric potential. The two cases investigated here, a copper shell and submersion in liquid argon, result in a similar reduction of the potential inside the detector underneath the passivated ring. However, the effect is stronger for the detector submerged in argon. The effect outside the detector is even stronger than inside, especially for the copper shell. However, the effect is particularly important for the case of submersion in liquid argon, where free radioactive ions can be attracted to the detector.  Alpha and beta decays of contaminants on the passivated surfaces of germanium detectors are known to be critical sources of background for rare event searches. Due to the special field conditions underneath these surfaces such events can look like 0 events. Thus, the knowledge on possible changes in the electric field and in the behavior of the detectors due to the surrounding medium, surrounding support structures and cables is important for the next generation of such experiments, which rely on a substantial reduction of the background level compared to what has been achieved so far.

Charge drift and induced signals
The drift velocity vectors are calculated separately for electrons and holes for each grid point using the electric field and the respective electron or hole drift velocity model. The default drift velocity models implemented in SolidStateDetectors.jl are based on the AGATA Detector Library [20,28] and are described in some detail in Appendix A.
For each energy deposition, the drift paths are calculated for point-like electron and hole clouds, using drift velocities, v /ℎ (r /ℎ ), as interpolated linearly from the neighboring grid points. The drift vector for each step is calculated as Δr /ℎ = v /ℎ (r /ℎ ) Δ , where the step length, Δ , can be chosen by the user and is usually in the range of 1-10 ns. The drift of a charge ends when it reaches an electrode or a maximum number of steps have occurred. If |Δr /ℎ | is below a certain limit the drift can also end inside the detector. For steps leading beyond the detector volume, the intersection between the drift vector and the detector edge is determined. If the crossing point does not belong to an electrode, the drift vector for this step is reduced to its component parallel to the surface. The charge drift ends when the crossing point belongs to an electrode.
Charges can get trapped inside the bulk of a germanium detector or, more commonly, underneath a passivated surface. The former is predicted by SolidStateDetectors.jl for volumes where the electric field is close to zero. The latter can, in general, not be predicted. There are several known effects like a surface charge-up or a surface channel, but there is no complete model for them. The user can, however, define virtual volumes, in which the drift model can be modified to simulate such effects. This is usually done after some effects have shown up in the data, which need to be understood. In many cases, a reduction of drift speed along the surface provides a first insight into what is happening.
The electric signals are the sum of the charges induced in the electrodes of the detector by all drifting electrons and holes. The Shockley-Ramo theorem [29][30][31] is used to calculate the time development of the induced charge, ( ), in each electrode from r /ℎ ( ) and Φ (r): where 0 is the absolute electric charge in the electron and hole clouds after separation of the charge-carrier pairs. Figure 5 shows the simulated drift trajectories and the signals induced in all electrodes for an energy deposition in the bulk of the n-type segmented BEGe detector. The separate contributions from electrons and holes show that the signal of the core electrode is dominated by the electron  drift in this type of detector due to the strong gradient of the weighting potential around the small core electrode, see Fig. 2c).
A positive collected charge is observed for segment 4, which is the collecting segment for this event. The other segments show mirror pulses that return to zero after the drift is completed. The amplitudes of the mirror pulses in segments 1 and 2 indicate that the position of the energy deposition is closer to segment 1 than segment 2. The possibility to obtain 3D position information from fits to the events in simulated pulse-shape libraries including all mirror pulses will be investigated in the future.
The parameters which are used to calculate the drift velocities, see Eq. (A.1) in Appendix A, are temperature dependent. The default values used in SolidStateDetectors.jl were determined for the reference temperature of 78 K [20,28]. However, germanium detectors are often operated at other temperatures. The effect of the temperature on the drift velocities can be taken into account by applying user defined model-dependent functions as described in Appendix B.

Comparison of simulation to data
Data from the n-type segmented BEGe detector, see Fig. 1, were used to validate the simulation. The detector was mounted with the core electrode on top in an electrically cooled aluminum vacuum-cryostat with a continuously monitored temperature control system [32]. The detector was irradiated in three measurement campaigns with 133 Ba, 137 Cs and 228 Th.
The signals from the core and segment electrodes were amplified with charge-sensitive preamplifiers and digitized with a sampling rate of 250 MHz. The core signal was inverted. The recorded pulses were baseline subtracted and corrected for their preamplifier channel-specific decay. Signal amplitudes were derived using a fixed-size window filter. Linear cross-talk between the core and the segments as well as between segments was corrected for in an automated energy-calibration procedure using calibration parameters obtained from single-segment events [24,33,34]. The calibration was performed for each data set individually. The data were not corrected for differential cross-talk. This kind of cross talk does not affect the measured energy but changes the shape of the affected pulses. A model based on the first derivative of the charge pulses, the currents, was developed [32] and applied to the simulation before a comparison to data .
Each experimental setup, including the respective collimator and source encapsulation, were implemented in G 4. The hits in the detector as simulated with G 4 were clustered and the induced charges were drifted using SolidStateDetectors.jl. The simulated pulses for the core and all segments were convolved with the response function of the respective preamplifier . The response functions were measured by injecting test pulses with a rise-time of a few nanoseconds as input to the preamplifiers and recording the resulting pulses. The decay of the preamplifier pulses was included in the response function and was varied within uncertainties to reproduce the variations observed in the data. Furthermore, the observed baseline noise was added to the simulated pulses. The resulting pulses were used for comparisons to calibrated and cross-talk corrected data.
The complete cross-talk model is that the vector of measured channel pulses, ( ), at each time, , is a linear combinations of all true pulses, ( ) and their first derivatives, ( ) = * ( ) + * ( )/ , where is the linear and the differential cross-talk matrix. Charges always drift to the core or a segment electrode; the segment boundaries do not cause any deficit in charge collection. The simulation is shown for charge injection at the average depth of 2.5 mm and a full G 4 simulation. The top right insert provides the segment numbering and the position of the interaction point on the 100 axis. Differential cross-talk was applied to the simulation [32]. The scales of the y axes for the mirror pulses were chosen such that remaining small differences between simulation and data become visible.

Surface events
Photons from the 81 keV line of 133 Ba were used to study events close to the detector surface. The data were collected by placing a collimated 133 Ba source on the side of the detector, which was kept at the reference temperature of 78 K. The pulses from single-segment events in the (81 ± 2) keV energy interval were averaged to form so-called super-pulses in order to reduce the effect of noise and to average out the spatial distribution of events due to the size of the beam spot. Figure 6 shows a comparison of simulated to measured super-pulses for a measurement where the middle of the mantle, = 20.2 mm, was irradiated at = 47.2 • . The super-pulses from simulation describe the data reasonably well. The core pulse is very well simulated. The simulated pulse of the collecting segment shows a slow-down at the end which is not observed in data and which is reflected in the mirror pulses as simulated for the non-collecting segments. The average penetration depth of 81 keV photons in germanium is approximately 2.5 mm such that the energy was deposited close to the surface and the holes were collected quickly. The electrons drift inwards along the 100 axis and upwards along the 001 axis. As the drift was dominated by electrons, the results shown in Fig. 6 validate the simulation procedure and the drift velocity model for electrons on the relevant axes at a temperature of 78 K.
The simulation is shown for the case where charges were injected at the nominal ( , ) position at a depth of 2.5 mm and the case where the charges induced by the individual hits as provided by G 4 were drifted separately. The isotope 133 Ba also features four gamma lines at higher energies (276, 303, 356 and 384 keV). These gammas penetrate deeper into the detector and can create background events with a deposited energy of around 81 keV through Compton scattering. The G 4 simulation included these gammas. As for the data, such events were included in the super-pulse formation. Nevertheless the super-pulses from charge injection and G 4 simulation are basically indistinguishable. This confirms that the signal to background ratio is high enough and that the detailed shapes of the initial charge clouds are effectively averaged out by forming superpulses for low energies. The remaining small differences between data and simulation observed for the segments could be due to the lack of knowledge about the exact impurity density distribution in the detector . They could also be a hint that the parameters of the electron drift model are slightly off. Such effects will be investigated in the future. Differential cross-talk was applied to the simulation and is unlikely to be a major source of remaining differences, see next section.

Bulk events
Data taken with a Compton scanner [32] using a collimated 137 Cs source mounted above the detector were used to compare simulated to measured pulses from energy depositions in the bulk of the detector. The Compton-scattered photons from the 662 keV line were detected by a cadmium zinc telluride (CZT) pixel detector with a position resolution of about 0.5 mm and an energy resolution of about 1% at the relevant energies. Events were selected which had coincident energy depositions in the BEGe and CZT detectors with a total energy of (662 ± 5) keV.
Events with two hits in the CZT detector were used to select energy deposits in specific regions of the BEGe detector within ( ± 1) mm, using a technique based on the reconstruction of Compton The impurity density distribution as provided by the manufacturer was multiplied with 0.9 to adjust the overall pulse length along the 100 axis. This is well within the uncertainty on the provided impurity density. No radial profile was implemented.
cones [32]. The and coordinates were taken from the source position. Super-pulses were formed by an energy-weighted average of the pulses from the selected events after normalization. Additional pulses from events with one hit in the CZT detector, which were compatible with the previous super-pulses, were added to the averaging for the final super-pulses of the chosen volume to further reduce the effect of noise in the data. Figure 7 shows the simulated super-pulses for an example location at = 25.5 mm compared to the data. The simulation agrees reasonably well with the data. This validates the simulation in the bulk of the detector, where the hole drift also plays an important role in the pulse formation. The data is corrected for linear cross-talk. The simulation is shown with and without the application of differential cross-talk [32]. Differential cross-talk is significant here because the time-structure of the pulses for the core and the collecting segment are sufficiently different. It mostly affects the segment pulses while the core pulse is basically not changed. The application of differential cross-talk improves the description of the data by the simulation, especially the description of the mirror pulses. The residual differences between the simulated and measured pulses of the collecting segment might be due to the lack of knowledge on the exact impurity distribution, which is an input to the simulation. They could also be an indication that the input parameters to the hole and the electron drift models need to be adjusted. The relative mobilities of holes and electrons are not very well known for germanium. Another effect to be investigated is diffusion and self-repulsion of the charge cloud. This will be implemented in SolidStateDetectors.jl in the near future.
Overall, the description of the data by the simulation is quite good. A quantitative statement can be obtained using pulse-shape analysis.

Pulse-shape analysis
One of the most common pulse-shape analysis techniques for point-contact detectors is based on the so-called / parameter, where is the maximum of the differentiated core pulse (i.e. current) and is the energy of the event. The performance of the simulation was tested using data, for which an uncollimated 228 Th source illuminated the detector from the side. It was placed 20 mm away from the surface at = 20.2 mm at the center of segment 1. Core pulses from a measurement lasting about 9 hours were used to measure the / distribution in selected energy regions. The simulated core pulses of events produced with G 4 for the experimental setup were smoothed and numerically differentiated exactly like the pulses recorded for data.
In Fig. 8, a comparison of the simulated / distributions to data and is shown for events in the double escape peak (DEP), single escape peak (SEP) and full energy peak (FEP) from 2615 keV photons from the decay of 208 Tl. The position of the peak in the / distribution derived from the 208 Tl DEP was scaled to 1 for data and simulation. The scaling factor for the data was only 0.5% larger than for the simulation. The simulated distributions describe the data quite well. However, the simulated distribution of / for the FEP shows events which have / larger than 1. Such events do not occur in data. They originate from large energy depositions directly underneath the core contact. This will be studied further.
The DEP from 208 Tl was used to determine a cut in / , which separates single-site from multi-site events. The FEP from 212 Bi at 1621 keV was used to test the performance of this cut. The spectra observed in data and simulation are shown in Fig. 9 for the relevant energy range together with the results of fits to the peaks, which were used to obtain survival fractions. When  The results of the study on / validate the simulation and encourage current and future 0 decay experiments using HPGe detectors to evaluate the signal efficiency of their pulse-shape discrimination with SolidStateDetectors.jl.

Developments
This paper represents the status of SolidStateDetectors.jl release v0.5. Substantial further developments are being incorporated in release v1.0, which is planned for the fall of 2021. The upcoming features can be classified as related to physics processes, tools to present the results or technical upgrades to further improve the ease of usage and the execution speed.
The implementation of the detector physics will be improved: • In release v0.5, the drift of individual point charges are simulated and combined by superposition. The interactions between them are ignored. Release v1.0 will have the option to drift charge clouds including the effects of self-repulsion and diffusion.
• Already in release v0.5, the electric potential and field as well as the drift velocity fields can be calculated for undepleted regions of the detector. In release v1.0, the weighting potentials will also be calculable for undepleted detectors.
• In release v0.5, the physics close to detector surfaces including the slow-down of charge carriers and trapping can be implemented by the user using virtual volumes. In release v1.0, a model will be implemented, where a charge cloud can be split and charges which come : Energy spectra from events from 228 Th a) data and b) simulation before and after the application of an / cut providing a survival rate of ≈ 90% for the 208 Tl DEP. Shown is the energy range around the 208 Tl DEP and the 212 Bi FEP. Survival rates were determined with fits performed with BAT.jl [35]. The confidence intervals 1 , 2 and 3 are shown for all fits. too close to a surface will be stochastically trapped. The user will have to provide input parameters like probabilities and what "too close" means in a configuration file. As these parameters are not really known the user might want to develop fits of predictions to data in order to determine them.
A number of tools in release v1.0 will help to better understand the detectors under study by providing information derived from the simulation: • Prediction of the full-depletion voltage; • Pulse-lengths for user defined pulse intervals, e.g. from 5-95% amplitude level; • Isochrones, i.e. surfaces from which charge drifts result in equal pulse-lengths; • Improved visualization capabilities.
The technical improvements encompass: • Improved constructive solid geometry (CSG); • Rotation, scaling and stretching of CSG primitives; • Definition of surfaces; • Distances to surfaces; • Support for fast calculation of fields on GPUs; • Usage of multi-threading beyond field calculations; • Extended documentation.

Summary
The open-source package SolidStateDetectors.jl has been introduced as a tool to simulate semiconductor detectors. While special emphasis is given to germanium detectors, the package can also be adjusted to simulate the response of silicon, or any other kind of semiconductor detector.
The package was programmed in the language ulıa with special care to ensure fast execution. Constructive solid geometry is used to define a detector and its surroundings. The package can handle three-dimensional calculations as required for segmented detectors as well as the usage of symmetries for simpler devices.
An n-type segmented point-contact detector was used to demonstrate the capabilities of SolidStateDetectors.jl and to successfully verify the simulation by comparing predictions to data. For rare event searches, accurate pulse-shape simulations can help to verify background identification techniques and provide guidance during detector manufacturing. A number of extensions for SolidStateDetectors.jl are planned for the near future. As an open-source package, SolidStateDetectors.jl is also open for improvements and extensions from the community.

Appendices A Drift velocity models
The drift velocity models were implemented for germanium with standard parameters applicable for standard detectors [32]. These models should, however, be directly applicable for silicon for the 100 and 111 axes, if the input parameters are adjusted. Further adjustments needed are provided as footnotes.
The drift of the charge carriers is governed by the electric field and the axes-dependent mobility tensors. The mobility tensors for electrons and holes are not fully known for germanium, but different parametrizations exist to describe the experimental results. As default, the drift velocity models from the AGATA Detector Library (ADL) [20,28] have been implemented in SolidStateDetectors.jl.
Germanium has a face-centered cubic crystal structure. When the electric field is aligned with one of the three principal crystallographic axes, 100 , 110 or 111 , the drift velocity vector is parallel to this axis. The longitudinal velocity, , describes the drift along these axes and is parametrized in the following way: where E is the electric field strength, while 0 , E 0 , and are parameters which were obtained by fits to measurements. The values of these parameters for electrons and holes for the 100 and 111 axes are set in configuration files. This allows for their optimization for a given detector. From the drift velocities along these two axes, the drift velocity of electrons and holes in any direction can be derived. For the calculations presented in this paper, default parameters [21,28] were used. In the following, the coordinate system is chosen with the axis, i.e. the symmetry axis of cylindrical detectors, aligned with the 001 axis.
A default to calculate the general electron drift velocity [36] is implemented. The conduction band in a germanium crystal reaches its minimal potential around the four equivalent 111 axes, called valleys . The electron drift velocity vector, v , can be written as where is the tensor of inverted effective masses for the th valley, / is the fraction of carriers in the th valley, E 0 is the normalized electric field vector and A (E) is a function of the magnitude of the electric field. The sum provides the direction of the electron drift. , where the rotation matrices are = (arccos √︁ 2/3) ( 110 + /2). The angle between the 110 and the axis, 110 , is set by the user for a given setup. The values = 1.64 and = 0.0819 [37] are given in units of the electron rest mass .
The deviation from an equal population of the valleys, / = 1/4 in germanium , depends on the electric field as [38] the components of the hole drift velocity vectors with Eq. (A.5). They are rotated from the local coordinates to the global coordinates as The predictions of this model have been verified [28] with experimental data [39]. There are indications that the description of mobilities may not be adequate for high-precision simulation of specific large volume detectors. If at all possible, the user is advised to adjust the mobility parameters for the detector design under study. Several functions , /ℎ ( ) are predefined. A specific function can be selected and its parameters can be defined in a dedicated configuration file. The scaling is defined for the two basic axes, 100 and 111 , and is propagated to all directions by the procedure described in Appendix A. The default is , /ℎ ( ) = 1 for all directions and electrons and holes, i.e. the default values are used for all temperatures.

B Temperature dependence of drift velocities
The most commonly used model for the dependence only considers the scattering of the electrons and holes off phonons in the crystal. It suggests a form of , /ℎ ( ) ∼ −3/2 [40]. This is qualitatively supported by some measurements for electrons where an , ( ) ∼ −1.6 behavior was observed. However, for holes, an ,ℎ ( ) ∼ −2.3 temperature dependence was observed [41,42]. Other simple analytical models exist for higher temperatures which also include the temperature dependence of the saturation velocity in high fields [43]. For some configuration, however, a Boltzmann-like temperature dependence of the drift time, −1 , ( ) = 0 + 1 − 2 / , was observed [44]; these data were incompatible with a power law. Any kind of functions , /ℎ ( ) can be introduced with the recommendation that , /ℎ ( 0 ) = 1 as long as the default parameters [20,28] for the drift models are used. The effect of the temperature dependence of the drift velocities can be neglected for many qualitative studies. However, the temperature dependence has to be taken into account for detailed comparisons of simulated pulses to data.