HPGe detector field calculation methods demonstrated with an educational program, GeFiCa

A review of tools and methods to calculate electrostatic potentials and fields inside high-purity germanium detectors in various configurations is given. The methods are illustrated concretely with a new educational program named GeFiCa - Germanium detector Field Calculator. Demonstrated in GeFiCa are generic numerical calculations based on the successive over-relaxation method as well as analytic ones whenever simplification is possible due to highly symmetric detector geometries. GeFiCa is written in C++, and provided as an extension to the CERN ROOT libraries widely used in the particle physics community. Calculation codes for individual detectors, provided as ROOT macros and python scripts, are distributed along with the GeFiCa core library, serving as both examples showing the usage of GeFiCa and starting points for customized calculations. They can be run without compilation in a ROOT interactive session or directly from a Linux shell. The numerical results are saved in a ROOT tree, making full use of the I/O optimization and plotting functionalities in ROOT. The speed and precision of the calculation are comparable to other commonly used packages, which qualifies GeFiCa as a scientific research tool. However, the main focus of GeFiCa is to clearly explain and demonstrate the analytic and numeric methods to solve Poisson's equation, practical coding considerations and visualization methods, with intensive documentation and example macros. It serves as a one-stop resource for people who want to understand the operating mechanism of such a package under the hood.


Introduction
The calculation of electrostatic potentials and fields in a high-purity germanium (HPGe) detector is the initial step in a full pulse-shape simulation [1][2][3][4][5] process. It is also used to guide the design of novel detector geometries to avoid unreasonably high depletion voltages or hidden undepleted regions, which may occur when the size of a detector is enlarged [6]. The design of read-out electronics can benefit from it as well since the capacitance of a detector, a determination factor of the electronics noise, can be calculated from the energy stored in the electric field in the detector. It is widely used in HPGe detector based neutrinoless double beta (0νββ) decay experiments, such as GERDA [7] and MJD [8], dark matter experiments, such as Co-GeNT [9], Texono [10] and CDEX [11], and gamma-ray tracking detectors to study structures of atomic nuclei, such as AGATA [12] and GRETA [13], etc.
SIMION is a commercial software package primarily used to simulate the transportation of charged particles in static or low-frequency RF fields. According to its documentation [16], it uses the finite-element method to calculate 2D and 3D fields with up to almost 20 billion grid points, given enough RAM. Its power in static field calculation is overkill for HPGe detectors while lacking some important features that are required for HPGe detector applications, such as the calculation of depletion voltage, region and detector capacitance, etc. This is understandable given that the main application of SIMION is not HPGe detector field calculation.
ANSYS Maxwell [17] is a more popular electromagnetic field simulation software compared to SIMION. It is for the design and analysis of electric motors, actuators, sensors, transformers and other electromagnetic and electromechanical devices. It uses automatic adaptive meshing techniques to achieve user-specified accuracy without detailed instruction from a user. As its main application is not HPGe detector field calculation, it has the same advantages and disadvantages as SIMION, but is more expensive than SIMION.
A common pitfall of all general-purpose commercial software is that one has to pay for extra features that are not needed in the HPGe field calculation, while still missing out some basic features that are needed.
FEniCS [19], on the other hand, is a free-to-use, open-source program developed by a global community of scientists and software developers, and is just as sophisticated as SIMION and Maxwell. Using efficient finite-element codes, its main purpose is to solve partial differential equations, including Poisson's equation, which is needed in HPGe field calculations. As versatile as it is, FEniCS demands effort to adapt it to a specific application, such as calculating fields in HPGe detectors. From this point of view, FEniCS has the same drawback as commercial packages, that is, it is overkill for HPGe field calculation, but lacks basic features that are specific for HPGe application. Nonetheless, there is ongoing effort within the MJD collaboration to adapt it for HPGe detectors.
On the contrary, MaGe [18] and siggen [1,14] are dedicated software for HPGe signal formation simulation. They are not as versatile and sophisticated as the previously mentioned packages, but are sufficient for the HPGe application. Initially, MaGe was jointly developed by the Majorana [8] and GERDA collaborations mainly as a GEANT4 [22][23][24] based Monte Carlo simulation package. It was extended later on to include a full pulse-shape simulation chain using GEANT4 simulation results as input [5,25,26]. It can be used for the simulation of both segmented [27] and point-contact [28] detectors. The major drawback of MaGe is that it is only available for the GERDA or Majorana collaborators.
Siggen [1,14] is mainly developed by David Radford for MJD pulse-shape simulation. It is open-source and free to use. A stand-alone portion of siggen, called fieldgen, is dedicated to the calculation of fields and potentials of point-contact detectors in two dimensional cylindrical coordinates. It cannot be used for segmented detectors. The program is written in c, but the configuration file is in plain ASCII with straightforward syntax for a user to easily specify detailed dimensions of a detector, such as the size of small electronic contact, or the width of a groove to reduce surface leakage current. Fieldgen can also be used to calculate the capacitance of a detector, the full depletion voltage, and the depletion region in case that a detector is not fully depleted. Those functions are not available in the packages mentioned previously. Fieldgen utilizes the successive overrelaxation method (SOR) to first calculate the potential in a coarse grid with a typical distance of 1 mm between two grid points. The result of this coarse calculation is then used as the input of a more precise calculation in a finer grid with a typical distance of 0.1 mm between two grid points. Using this simple approach in place of automatic adaptive meshing techniques used in some of the other packages makes fieldgen both fast and accurate enough for its dedicated application.
SSD [20] is mainly developed by the GeDet group at the Max-Planck-Institut für Physik for LEGEND [21]. It is capable of not only the calculation of electric fields but also the simulation of electronic signals. In its field calculation part, it contains functions to deal with common detector configurations, such as point-contact and segmented ones. It features adaptive grid sizes, which improves both the calculation speed and accuracy. It is written in Julia [29], a relatively new, high-level, general-purpose programming language designed to address the needs of high-performance numerical analysis. It is possible to enable multi-threading in SSD, which takes the full advantage of modern computer hardware. Compared to MaGe and fieldgen, SSD has an attractive feature to calculate the field outside of a detector taking into account the influence of the detector holding structure nearby.
Overall, fieldgen seems to be the maturest at this moment for users interested in HPGe field calculation as long as their detector geometry is similar to that of point-contact ones 1 . However, the lack of detailed documentation makes it hard for a developer to modify the code of fieldgen for other geometries or to add new features. This is one of the reasons why many research groups write their own code for HPGe detector field calculation instead of using the mentioned major players. An obvious advantage of home brewed code is that it is well understood and easy to tune if needed. The second advantage is that writing their own code instead of using existing ones deepens the understanding of junior researchers on HPGe detector working principles and numerical calculation techniques. Drawbacks of this approach include the limited functionality, the lack of verification and the waste of time in reinventing the wheel.
GeFiCa is aimed at clear explanation and demonstration of the analytic and numeric methods to solve Poisson's equation, practical coding considerations and visualization methods. It does so by providing intensive documentation and example macros, and serves as a one-stop resource for people who want to understand the operating mechanism of such a package under the hood. None of the tools mentioned above fits all applications. Home brewed codes built on top of some existing tools may be the best choice for education and specific applications, as long as the drawbacks mentioned previously can be effectively overcome through the demonstration provided in GeFiCa.

Space charges
HPGe crystals come in two types. As shown in Fig. 1, if the trace impurity atoms in a crystal provide freemoving electrons (phosphorus, for example), the crystal is of n-type, and if the atoms provide free-moving holes (boron, for example), the crystal is of p-type. In both literature and popular science articles, these freemoving charge carriers are often preceded with adjectives like "extra" or "excess", which may lead to a false impression that an n-type crystal has "extra" electrons donated by donor impurity atoms and is hence negatively charged, or that a p-type crystal has "extra" holes (vacancies in covalent bonds) due to acceptor impurity atoms and is positively charged. These free-moving charges can be regarded as "extra" since they are not used in forming covalent bonds between atoms, which is the fundamental reason why they are free. But they are not "extra" charges that break the balance of the numbers of protons and electrons in a crystal. Actually, no matter which type it is, a crystal is electrically neutral as a whole because the number of protons are the same as the number of electrons in both impurity and Ge atoms. As trivial as it sounds, this fact is worthy of emphasizing, especially for one to understand the sign of space charges to be mentioned in the following paragraph. When the bias voltage applied to a crystal is high enough, all free-moving charge carriers can be swept out of the bulk of the crystal. The crystal is said to be depleted of free charge carriers. In an n-type crystal, it is the free-moving electrons that are swept out. Consequently, the trace impurity atoms are positively ionized. Since the ions are fixed in their locations in the crystal, they cannot be swept out by the external electric field, and are hence called "space charges". In a p-type crystal, however, the space charges are negative, since it is the free-moving holes that are swept out. It is quite counter intuitive for one to realize the fact that a depleted n-type crystal is actually positively charged and a p-type negatively charged. Space charges create an electric field in addition to the one that is created by the bias voltage. The total electric field inside a depleted crystal is a linear combination of these two.
The space charge density distribution, normally denoted as ρ, can be quite complicated due to the nature of HPGe single crystal growth process [30,31]. It is normally characterized in the following way. First, a few wafers are cut from various axial positions in a HPGe single-crystal boule pulled using the Czochralski method, typically, one from the shoulder and one from the tail of the boule. Second, small samples are cut from individual wafers along their radius. Net impurity concentrations of these samples, N A − N D , are then measured using Hall-effect, where N A is the acceptor concentration, N D donor concentration, with a unit of cm −3 . Since there is a relationship between ρ and N A − N D as explained in the previous paragraph: where e = 1.6 × 10 −19 C is the elementary charge, both the vertical (axial) and radial distributions of ρ can be investigated this way. A typical vertical net impurity concentration profile of a HPGe single-crystal boule is shown in Fig. 3 taken from [30] with a vertical double-dotted dashed line added to clearly indicate the p-type and n-type regions. The dashed line indicates the contribution from a typical p-type impurity element, Al. The dotted dashed curve indicates the contribution from another typical ptype impurity element, B. The dotted curve shows the contribution from a typical n-type impurity element, P. The solid curve broken around 80% of the boule is the overall net impurity concentration. The crystal is of p-type from 0 to 80% of its length, and changes to n-type after that. The curve is approximately flat from 20% to 40%, which is a typical portion of the boule to be harvested for detector fabrication. A typical radial net impurity concentration profile of a HPGe single crystal is shown in Fig. 4 taken from [31]. It is basically flat from 0 to a certain radius, but increases dramatically close to the skin of the crystal. Sometimes, a crystal may even change its type from its center to its outer radius, as mentioned in [30]. The skin of a boule may be removed so that the central part used for detector fabrication has a relatively constant impurity distribution. Given those experimental evidences, the space charge density in general has to be expressed as a function of location, i.e., ρ(x), where x is a vector indicating the location of interest. Since the measurement of impurity is destructive for the raw material, the real impurity distribution in a crystal used for detector fabrication is usually unknown. Normally, only the average impurities close to the top and bottom of the cut portion of the crystal are known. The impurity distribution in between is regarded as a constant or approximated by a first-order polynomial determined by the average top and bottom impurities. If the right portion of a crystal (20-40% of the black line in Fig. 3, for example) is harvested for detector fabrication, this is normally an acceptable approximation. However, one has to keep in mind that our knowledge of the real impurity distribution is incomplete, and our approximation may have sizable uncertainties.

Poisson's Equation
The existence of space charges complicates the calculation of the electrostatic potential in a HPGe detector. Without space charges, the potential can be calculated by solving Laplace's equation, where V (x) is the potential to be determined. With space charges, however, the potential must be calculated by solving Poisson's equation, which takes into account the space charge distribution in the bulk of a detector: where = 0 r with 0 ≈ 8.854 × 10 −12 F/m being the permittivity in vacuum, and r ≈ 16.0 being the relative permittivity (or dielectric constant) of Ge. Both differential equations have an infinite amount of solutions characterized by a few undetermined constants. These constants can be fixed by boundary conditions, which refers to the voltage values on detector electrodes. The relationship between these two equations can be understood better when we consider two different boundary condition setups: first, potentials of electrodes of a detector are set based on the bias voltage applied to the detector, second, they are all set to zero. If Laplace's equation is solved with the first setup, its solution is a potential field caused by the bias only. If Poisson's equation is solved with the second setup, its solution is simply the potential caused by the space charges only. The potential in a detector is a linear combination of these two solutions. We can also solve Poisson's equation with the first set of boundary conditions, which directly results in the combined potential.
In addition to the potential, we are also interested in the electric field distribution in a detector. The electric field vector E can be then determined with the equation These equations are rather abstract. A concrete expression can be obtained in a specific coordinate system. For example, in Cartesian coordinates, Poisson's equation reads, where, x, y, z are the three Cartesian coordinates. The expression of Poisson's equation in spherical and cylindrical coordinates are listed in Appendix A.

Planar Detectors
As mentioned in section 2, in general, the space charge density ρ is a function of x, and there is no analytic solution for three dimensional Poisson's equation with a complicated ρ(x). However, in certain highly symmetric detector configurations, Poisson's equation can be significantly simplified and its analytic solution can be obtained. For example, at the center of a large but thin planar HPGe detector, the electric potential can be regarded as only varying along the thickness of the detector, x, and ρ can be regarded as a constant. Eq. 5 can then be simplified to Its analytic solution reads, where C 1 and C 2 are constants which can be determined using two boundary conditions, i.e., the voltages of two planar detector electrodes. The electric field reads,   6 shows the voltage versus the thickness of a planar detector, assuming a thickness of 1 cm and a voltage of 2,000 V applied to its top electrode. The net impurity concentration corresponding to each curve in the figure is listed in the legend. When ρ = 0, Eq. 7 becomes V = C 2 x + C 1 , which is simply a straight line between [0,0] and [1 cm, 2,000 V]. The higher an impurity concentration, the more a curve is bent up or down depending on the type of the impurity. Since the slope of curves in Fig. 6 is proportional to the magnitude of the electric field, as shown in Eq. 8, the bending of the curves shows how space charges modify the overall electric field in a detector. This is demonstrated explicitly in Fig. 7, where the y-axis changes to electrical field in the unit of V/cm. A small change of the impurity concentration may result in large deviation of the overall field from the constant external field. When the impurity concentration is high enough, the electric field close to the electrodes of the detector can be as low as zero. Such low field regions are where severe charge trapping may happen, which deteriorates the energy resolution of the detector, hence are not desirable. Obvious solutions of such a problem include applying a voltage significantly higher than the depletion voltage, reducing the thickness of the detector, or growing purer crystals. The first solution is dangerous, the second is undesirable and the last is difficult. A less obvious alternative is to switch to a different geometric configuration of the detector.  How the geometry of a detector can help solve this problem can be clearly demonstrated using the ana-lytic solution of Poisson's equation in cylindrical coordinates. The electric potential far away from its two end surfaces of a true-coaxial HPGe detector can be regarded as varying only with r. If one further assumes that the space charge density ρ is a constant, Poisson's equation in cylindrical coordinates can be simplified to

Coaxial Detectors
Its analytic solution reads, where C 1 and C 2 are constants, which can be determined using boundary conditions, that is, the locations and voltages of the two electrodes of a detector. The electric field is then Fig. 9 Voltage versus radius of a true coaxial detector. Fig. 9 shows the voltage versus the radius of a truecoaxial detector with an inner radius of 0.25 cm, an outer radius of 1 cm and a voltage of 2,000 V applied to its inner electrode. The net impurity concentration corresponding to each curve in the figure is listed in the legend. Compared to Fig. 6 for a planar detector, the curve corresponding to the zero impurity is not a straight line anymore. Instead, it bends downward, reflecting the fact that the electrical field close to the inner radius is stronger, which can be seen in Fig. 10 as well. Now, p-type impurities bend the curves further down, while n-type impurities bend them upwards, effectively flatten the electrical field distributions along the radius, as shown explicitly in Fig. 10. Given the right impurity concentration, the electrical field in a true-coaxial detector can be optimized to avoid charge trapping. To demonstrate this point more clearly, a curve corresponding to a high impurity concentration of −6×10 10 /cm 3 is added to Fig. 9 and 10, which is not in Fig. 6 and 7. The electric field distribution corresponding to this concentration is in between ∼2,000 V/cm and 4,000 V/cm, well above zero in the entire volume of the detector. One has to avoid falling into a false impression that coaxial detectors prefer n-type crystals to p-type ones. In reality, this preference can be easily flipped by flipping the bias polarity. It is better to say that the type of the crystal prefers a certain bias polarity. As a conclusion, coaxial detectors are much more tolerable for high impurity concentrations than planar ones.  In spherical coordinates with polar and azimuthal symmetries, the θ and φ terms can be dropped and the Poisson's equation can be simplify to

Hemispherical Detectors
Assuming constant ρ, its analytic solution reads, where C 1 and C 2 are constants that can be determined by boundary conditions. The electrical field A detector with such a configuration is more impurity tolerable than a coaxial one. However, it is not easy to bias the inner radius of a full sphere. The electrical field of a hemispherical detector can be approximated with the same solution, and it is possible to apply voltage to its inner radius. However, it is a significant challenge to machine a cylindrical single-crystal boule into such a shape. For this reason, no hemispheric HPGe detectors have been made so far, and most HPGe detectors take the cylindrical shape for convenience. If the inner radius of an imaginary hemispheric detector is small enough, say, about 1 mm, it can help illustrate some important properties of a point-contact detector, which can be imagined as a traditional coaxial detector with its central contact shrunk to a point. For example, certain amount of impurity is necessary to shape the electrical field in the detector so that it is not too strong close to the point-contact and not too weak far away from it. This is why the first point-contact detector is called a "shaped-field" one [28].

Depletion Voltage
Given fixed dimensions and impurity concentration of a crystal, we'd like to find out the voltage at which the crystal can be fully depleted, or the depletion voltage, V d . The method to solve this problem can be demonstrated using the analytic solution, Eq. 7, of the onedimensional Poisson's equation in Cartesian coordinates. The strategy can be applied to multi-dimensional configurations with minor modifications.
To keep our discussion as concrete as possible, let us assume an ideal planar detector with a thickness of d = 1 cm and a homogeneous impurity concentration of 4 × 10 10 /cm 3 (p-type) in its entire volume. Let's further assume that its bottom electrode is at x = 0 and grounded, i.e.
Applying this boundary condition to Eq. 7, we have C 1 = 0. If no bias is applied at the top electrode, that is, V (x = d) = 0, we can further get where ρ = −4 × 10 10 /cm 3 e is the space charge density. Eq. 7 can then be rearranged as, which is the green parabola shown in Fig. 12. However, this is not physically correct, since the whole crystal should be at V = 0 without any bias. The problem comes from our taken-as-granted assumption that ρ = −4 × 10 10 /cm 3 e, which is only true in depleted region.
In undepleted region, there should not be any space charge, that is, ρ = 0, which indeed guarantees V (x) = 0 in Eq. 16. Instead of regarding it as a mistake, there is a better way to interpret Eq. 16, that is, it is the contribution to the voltage from the "space charge alone", when the detector is fully depleted. Its value hence is not dependent on the bias voltage after fully depletion. The overall voltage should be the sum of this contribution and the voltage due to the external bias.
The external bias voltage distribution without any contribution from space charges is simply V (x) = C 2 x (Eq. 7 with ρ = 0, C 1 = 0), that is, a straight line, as the one labelled "bias alone" in Fig. 12. At a bias voltage of −1000 V, the sum of "space charge alone" and "bias alone" contributions gives the red curve in Fig. 12 that is below the −1000 V line in a wide region. This region can be regarded as a potential well where positive free charge carriers are trapped. In another word, the −1000 V bias is not enough to sweep all free charges out of the crystal. The curve is hence labelled "undepleted".
What we are interested in here is to identify differences between an undepleted case and a depleted one. In this concrete example, an "undepleted" curve has |V (x)| > |V d | at some x. If we change the crystal from p-type to n and keep other configurations unchanged, the "space charge alone" curve would bend downward. More general criteria hence would be, an "undepleted" curve has V (x) out of the range defined by boundary voltages, or dV / dx changing its sign, at some x.
Assuming a certain bias voltage and a constant ρ over the whole volume, if the final answer is "undepleted" according to the criteria identified previously, we have to start over again assuming a higher bias. Obviously, the detector will certainly be depleted given an extremely high bias. However, in reality, it is hard to deliver a very high voltage without micro (or even major) discharges along the high voltage cable. Normally, the operation voltage is ∼ 1000 V over the depletion voltage.
The difference between an over depleted curve and a just depleted one, in this concrete example, is that for the latter, but E(x = d) > 0 for the former case.
In general, we have to do a search in between 0 and a large bias voltage for the "just depleted" case, where the electric field E on one of the boundaries is exactly zero. The calculation for this analytic example is very fast. Special treatment has to be taken in multi-dimensional numerical calculations to avoid expensive computations (See Sec. 5.3).

Impurity Requirement
Given technical difficulties in delivering high voltages in a cryogenic environment, a low depletion voltage is always preferred. It is a common practice to figure out the maximal net impurity concentration a crystal with certain dimensions must have to be depleted at or under a given voltage. In our previous example, the depletion requirement is E(d) = 0, which allows us to calculate C 2 : Insert the calculated C 2 and C 1 = 0 back to Eq. 7, we get the maximal allowed space charge concentration ρ = 2V /d 2 , where V is the given voltage. Analytic solutions are not available for more complicated detector configurations. In that case, we need to make guesses on the impurity concentration or even profile, search for corresponding depletion voltages based on the method described in Sec. 5.3, and see if they go below the required voltage.

Numerical Calculation
Even though many detector design concepts can be demonstrated with analytic solutions of highly symmetric detector configurations, numerical calculations are necessary for more advanced configurations that cannot be simplified to lower dimensional problems.
The first step of numeric calculation is to establish a grid within the detector volume, which consists of many tightly spaced points, some right on boundaries, others inside. The field values of a grid point can be determined by those of its immediate neighboring points. Their relations are dictated by Poisson's Equation in its numeric forms. Starting with the known values of the points on boundaries, the value of each point can be uniquely determined.
Configuring a grid that ensures an efficient and accurate calculation is an art by itself. For the sake of clarity in our discussion without losing generality, let's at first consider a section of a one dimensional (1D) grid around a point at x, as shown in Fig. 13. Numerically, the second order derivative on the left side of Eq. 6 can be expressed as where dx ± are the distances from the point at x to the previous and the next points as shown in Fig. 13. It is possible to involve more points in the calculation, such as the previous previous or next next points, but the basic idea is the same. There are different ways to rearrange Eq. 6 based on Eq. 18, which lead to different methods to solve the problem. The two most common ones are the conjugate gradient method and the successive over-relaxation method.

Conjugate Gradient Method
The conjugate gradient method starts by moving all known terms, such as the boundary voltages and terms containing ρ(x), etc., to the right side of Eq. 6. Assuming n points in our 1D grid shown in Fig. 13, and V (x) = V i , ρ(x) = ρ i are the values at the i th point, Eq. 6 becomes, We have such an equation for n − 2 points, excluding the first and n − 1 one, since V 0 = 0 and V n−1 = the bias voltage are known and have to be included in K i 's. The n − 2 linear equations can be collectively written as where C is a n − 2 by n − 2 matrix, and V and K are vectors with n − 2 elements. C is sparse, with at most 3 non-zero elements in each row. It is also symmetric and positive definite. A standard way to solve such a linear equation system is the conjugate gradient algorithm [32], which boils down to minimizing a quadratic function of V with the form V T CV /2 − K T V . There is a ROOT [33] macro included in GeFiCa to demonstrate the method. It works well when n is below 100, but becomes painfully slow for a large n.

Successive Over-Relaxation Method
Another way to rearrange Eq. 6 is If ρ = 0 and dx − = dx + , it can be simplified to where V i is simply the average of its neighboring values. In both equations, V i is calculated given V i−1 and V i+1 . However, since V i−1 and V i+1 are also unknown (except for V 0 and V n−1 ), we need to start the calculation with some initial values. One choice would be V n−1 = the bias voltage, V bias , where the superscript (0) indicates that these are the initial values of grid points.
Given these initial values, we can use Eq. 20 or 21 to update V i . Use Eq. 21 in our calculation hereafter to simplify the demonstration, we have where j indexes the steps of updating. Since V (j) it pulls the value of its neighbor V n−2 up a bit after each updating, and V n−2 pulls up V n−3 , and so on and so forth. After many iterations, V i becomes very close to its true value, the difference between the values of V i in current and previous iteration becomes very small. We can use the following criterion to stop the iteration: < a small value, e.g.,10 −8 . (23) This is the so-called Successive Relaxation (SR) method.
To speed up the relaxation process, we can manually increase the difference of a value between two iterations by introducing a constant, 1 ≤ F R ≤ 2, in the following way: where, V , is the value updated by the original relaxation method, F R is called the relaxation factor. This is why the method is called Successive Over-Relaxation (SOR) method. The concept of SOR is depicted in Fig. 14, where the first a few steps of updating are shown for the last a few grid points in an idea planar detector without any impurity. A carefully chosen relaxation factor can reduce the total number of iteration significantly, which is discussed in detail in Sec. 8.1. The same method can be applied to multi-dimensional problems in various coordinate systems. The counterparts of Eq. 20 in those systems are summarized in Appendix B.

Depletion Voltage
The general method described in Sec. 4.4 applies to numeric calculations as well. We need to search for a V d that just depletes the detector with the following procedure: 1. Pick up a V min , say zero, and a V max , say 10 6 V. 2. Assume a bias voltage, V bias ∈ (V min , V max ).
One feature of the undepleted curve in Fig. 12 distinguishes it from the depleted ones; that is, the maximum or minimum of the potential is not on the boundaries of the detector. Inspired by this, the criterion of depletion in step 4 for numerical calculations can then be set as none of the grid points has a potential that is larger or smaller than the value of any of its neighboring points.
A potential drawback of the described method is that it may be time consuming if every new search needs to run an SOR. Fortunately, there is a way to avoid that. As demonstrated in Fig. 12, the total potential distribution, V i , is a sum of the distribution due to impurity alone, V ρ i , and the one due to bias alone, represents the potential distribution due to unit voltage, 1 V, the total potential distribution can be calculated as If V u i and V ρ i are calculated using SOR before the described iteration, step 3 can be replaced by Eq. 25 instead of another SOR.

Undepleted Region
When |V bias | < |V d |, some region of the detector is not depleted. Numerically, the undepleted region can be found by applying the following procedure to every grid point in the SOR process: 1. Calculate the potential of a grid point V i using potentials of its immediate neighboring points. 2. Find the maximal and minimal potentials V max and V min of the immediate neighboring points. 15 shows potential distributions of a planar detector with a |V bias | < |V d | after some chosen numbers of SOR iterations. One can see how the undepleted region grows larger near one of the electrodes. Another interesting thing to notice is that it does not take many iterations for the potential to become very close to its final values. Most iterations after that are used to improve the accuracy in a few percent level. As shown in Fig. 15, the undepleted region in a planar detector is adjacent to one of its electrodes. In the case of a point-contact detector, the undepleted region can stand alone somewhere in the center of the detector, away from any electrode. This is the so-called pinch-off effect, since the depleted region is "pinched off" from electrodes.
Since there is no electric field in the undepleted region to separate and drift electrons and holes generated by radiation interactions, the region is insensitive to radiation. Even if a pair of charge carriers are generated in the depleted region, one of them may drift to the undepleted region along the electric field and get stuck there instead of being collected by an electrode. It is hence worthwhile to reveal the existence of such an undepleted region through field calculation. Fig. 16 shows the electric field as a colored contour in logarithmic scale in a point-contact detector. The point contact is at the origin of the plot. The field value around the center of the detector is very close to zero, so the logarithm of them approaches negative infinity. They are color coded as white, which nicely visualizes the undepleted region that is pinched off from the point contact.
This phenomenon seriously limits the size of a pointcontact detector, since the electric field inside the detector becomes weaker when the size of the crystal becomes larger if the bias is not ramped up accordingly. To avoid this, one can bore a central hole from the opposite side of the point-contact, metallize its surface and keep it at the same bias as other surfaces. Such a detector is called a inverted-coaxial point-contact detector, or ICPC in short [34]. Fig. 17 shows the electric field distribution in an ICPC as a color coded contour in logarithmic scale. Other configurations of this calcu- lation are kept the same as the ones used to generate Fig. 16, including the crystal impurity level and the bias voltage. One can see clearly that the undepleted region is successfully eliminated from the center of the crystal.

Electric Field Lines
The thick black lines in Fig. 16 and 17 are estimated charge drift trajectories starting near the outer surface of the detectors. The procedure of the estimation can be illustrated in two dimensional Cartesian coordinates: 1. Linearly interpolate the electric field components E x , E y at a random starting point (x, y) using values at its four neighboring grid points. 2. Calculate the total electric field E = E 2 x + E 2 y at the same point. 3. Calculate the propagation of a positive unit charge along x and y: dx = µE x dt , dy = µE y dt, where µ takes a value of 40, 000 cm 2 /(volt·second), a number in between typical electron and hole drift mobilities in HPGe crystals [35][36][37], and dt takes a value of 10 ns. 4. dx, dy are then further modified using equations dx = dx × weight, dy = dy × weight, where weight = (5 volt/mm)/E, which is used to stretch dx , dy in a weak field and shrink them in a stronger one. 5. The new position of the positive unit charge is then updated to (x + dx , y + dy), which is saved in an object of the FieldLine class. 6. Repeat step 1 to 5 using the updated starting point coordinates until it moves out of the crystal or falls into an undepleted region.
Changing the positive unit charge to a negative one would let it propagate to the opposite direction.
Ignoring the influence of the germanium crystal structure, charge carriers drift roughly along electric field lines. The propagation path created this way can hence be regarded as a rough estimation of the field line.
It is interesting to see in Fig. 17 that the field lines merge in the middle of the detector and get collected at the point contact, just as streams flow down to a river in a valley (the blue-ish region in Fig. 17 if color-printed).

Boundaries in between Grid Points
Sometimes, the edges between the side and end surfaces of a cylindrical detector are tapered, shown as the small white triangular regions at the corners of the color contours in Fig. 16 and 17. A crystal boundary line hence can go in between grid points that are distributed along orthogonal lines, as shown in Fig. 18. Assuming a simple case, where the grid points are evenly spaced, the distances between them take fixed values, dx and dy. The distances of a regular grid point to its previous and next neighbors equal to each other: dx − = dx + = dx. In case of a point near the boundary, such as (i, j) shown in Fig. 18, it is more precise to replace dx − with dx − i when evaluating Eq. 18 along the x-axis, where dx − i is the distance to the boundary instead of the distance to the previous grid point as shown in Fig. 18. Similarly, dy + should be replaced by dy + i for a more precise evaluation of Eq. 18 along the y-axis. This effectively moves nearby points on the vacuum side right to be on the boundary, shown as the red dots in Fig. 18. Such an operation can only be done if variable step lengths are allowed for individual grid points.

Weighting Potential in Segmented Detectors
In addition to the real electric field, the so-called weighting potential is also of great interest, since it can be used to calculate the electric charges on an electrode induced by the drifting charge carriers inside a detector based on the Shockley-Ramo's Theorem [38][39][40]. It differs from a real potential in two ways. First, it is purely determined by its boundary conditions. The impurity concentration in a crystal should be regarded as zero in calculating the weighting potential. Second, the voltage on the interested electrode should be set to 1 volt, while the voltage on any other electrode should be set to zero in calculating the weighting potential. The weighting potential of an electrode is probably best demonstrated using a segmented HPGe detector. Fig. 19 shows the cross section of a detector segmented evenly in six along the azimuthal direction. The weighting potential of one of the segment electrode is overlaid as a colored contour in a logarithm scale. The white circle in the middle indicates the core electrode of this cylindrical detector. The colored contour does not quit reach the bottom boundary, simply because the potential there is too close to zero to be color coded in a logarithm scale.

Capacitance
The capacitance of a HPGe detector C d is of special interest due to at least two reasons. First, the electronic noise of a HPGe detector is proportional to the sum of C d and the capacitance of the feedback capacitor, C f , in the pre-amplifier circuit of the detector [39,41,42]. Second, C d decreases as the detector bias voltage ramps up. The reason becomes clear later in this section. This feature can be used to measure the depletion voltage, V d . It can also be used to check if a detector operates normally during the ramping up of its bias voltage. It is therefore an important task of a field calculation package to calculate C d given an arbitrary bias voltage. For an ideal one dimensional planar detector, where A is the area of an electrode of the detector and d is the thickness of the depleted region in the detector, which can be calculated as This relation can be derived from Eq. 7 with the boundary condition 15 and 17. C d hence is anti-proportional to √ V d , and decreases as V d increases, until d becomes the thickness of the plan detector. After that, C d stays at its minimum since d cannot increase anymore. The square data points in Fig. 20 are calculated using Eq. 26 and 27 given individual bias values.
The depletion depth d can also be determined numerically using the method described in the previous section. C d can be then calculated using Eq. 26. The results are the triangle data points in Fig. 20.
For a detector configuration as complex as a pointcontact one, there is no analytic solution for C d . The following numerical method is used in GeFiCa to calculate C d . It is based on the fact that the energy stored in a charged capacitor U is equal to the overall work done W to move a total amount of charges Q to the electrodes against the electric field E caused by Q stored in the capacitor: Given an arbitrary amount of charges q already stored in a capacitor, the work done to increase it by an infinitesimal amount dq is Integrating it on both sides yields The relation C d = Q/V bias is used in the last step of the derivation to replace Q, an unknown variable, with C d and V bias . On the other hand, since the electric field energy density is E 2 /2, U can be expressed as where dτ is the volume integration element. For a planar detector with a constant impurity, the integral can be solved analytically as: Replacing U and W in Eq. 28 with Eq. 32 and 30, we derive Eq. 26. The numerical version of Eq. 31 for an ideal planar detector in Cartesian coordinates is where i is the index of each grid point. Combining Eq. 33, 30 and 28, C d per unit area A can be calculated as This is implemented in function GeFiCa::X::GetC().
The results are shown as the circle data points in Fig. 20. The perfect agreement between all methods verifies two numerical calculations in GeFiCa: the finding of the undepleted region (or depleted region) and the calculation of C d given an arbitrary V bias . It is worth noting that the electric field E here is only due to Q accumulated on the detector electrodes. It is different from the actual field in a depleted detector which is the combination of the fields from both Q and the space charge in the crystal.
For a point-contact detector in Cylindrical coordinates, the numerical version of Eq. 31 is It is implemented in function RhoZ::GetC().

Interpolation Between Grid Points
A numeric calculation can only give the field values right at each grid point. Interpolations are needed to get the field values at a random point that may not coincide with any grid point. Fig. 21 shows the equations to linearly interpolate the potential value at the point of interest that falls in between four points in a 2D Cartesian grid using the known potential values on those four points, V 1 , V 2 , V 3 and V 4 , taking into account the distances between points, X, Y, x, y. This method does not work for unit grid squares across crystal boundaries that are neither in parallel with nor perpendicular to grid lines, since those boundary lines can separate a square into irregular shapes, the interpolation of which can be complicated. There are three ways that such a boundary line can go through a unit grid square as shown in Fig. 22. Potentials at the crossing point, V A and V B , are equal to the bias applied to that boundary. Most of the time the field outside of the crystal is not of interest. Within the crystal, the point of interest can fall into either a triangular or a rectangular region marked as T or R, separated by blue dotted lines in Fig. 22. If it falls into an R region, the interpolation method shown in Fig. 21 can be used. If it falls into a T region, the triangular interpolation shown in Fig. 23 can be used.
where (x 1 , y 1 ), (x 2 , y 2 ) and (x 3 , y 3 ) are the Cartesian coordinates of the grid points 1, 2 and 3 shown in Fig. 23, (x, y) are the Cartesian coordinates of the point of interest. In case of the vector field E, interpolations are done separately for individual components to get E x , E y at the point of interest. The total E is then calculated as E 2 x + E 2 y .

Coding Convention
The coding convention is similar to that of ROOT [43]. The following exceptions are used to increase the readability of the codes to the user: -Class names do not have prefix letters, such as T in ROOT. Instead, the name space GeFiCa is used to avoid name collision should GeFiCa be used together with other libraries. -Configurable member variables are made public to avoid trivial getters and setters. Their first letters are capitalized. Unlike private member variables, they do not have letter f prefixed.

Class Structure
As shown in Fig. 24, most of the GeFiCa classes belong to two categories: grid and detector. Those that are derived from class Grid are used to describe grid setups. Those that are derived from class Detector are used to describe detector configurations. The Grid class inherits a set of arrays from the Points class to describe variables associated with individual grid points, such as coordinates and field values. Names of its derived classes indicate the dimension and coordinates used to construct the grid. For example, X is used for a one dimensional grid in Cartesian coordinates, RhoZ is used for a two dimensional grid in cylindrical coordinates. The Detector class inherits impurity setup from the Crystal class. Its derived classes, such as PointContact, inherit from it the common detector setups, such as bias voltages. A grid class can get boundary conditions and the impurity distribution from a corresponding detector class through a virtual function interface defined in the Grid class: The data flow can be the other way around, that is, a detector class gets grid setups from a grid class. However, since it is the grid that the SOR process updates instead of the detector configuration, this is a less natural choice. With the current data flow direction, the SOR can then be performed by simply calling Another choice would be to combine the detector and grid classes. For example, instead of having both PointContact and RhoZ, we can create a single class called PointContactRhoZ. The advantage of this approaches is that there is no need to pass information from the latter to the former through some interface functions. The drawback is the lack of clarity, the same class object will be used for both detector configuration and grid operation. Considering the main purpose of GeFiCa is to demonstrate the logic, methods, and techniques for field calculation, we chose not to use this approach.
To its root, this is actually a question of to what extend we want to utilize the object-oriented (OO) coding style. Think about two extreme cases. First, we can write everything in a single main function. Second, we can create a class for each individual functionality, such as the impurity profile and the bias voltage. The first approach relies on careful documentation to clarify its internal logic. The second introduces many trivial interfaces to pass information between classes. A balanced approach in between is adopted for GeFiCa.

Data Structure and I/O
Minimally, two float numbers are needed for each point in a one dimensional grid with a fixed interval between its points: one for the spacial coordinate and another for the electric field potential. The number of points in a grid must be changeable according to the dimension of a detector and the precision of a calculation. This demands the use of arrays that can change in size to store variables of individual points. Even though a float number is precise enough to hold the final result of a numerical calculation, a double is preferred to preserve precision during iterations of a SOR process. A standard C++ vector<double> is used in GeFiCa for each variable to provide enough precision and flexibility.
Given a grid with variable step sizes as shown in Fig. 13, one more variable is needed for each point to store the distance to its next neighbor, dx + i . The distance to its previous neighbor is saved as dx + i−1 in its previous point. In GeFiCa, however, the distances to both the previous and next neighbors, dx − and dx + , are saved, since they may be different for some of the boundary points as detailed in Sec. 5.6. This certainly creates redundancy in storage for points away from boundaries. However, since output files of GeFiCa are saved in ROOT format instead of plain ASCII, such redundancy does not increase their sizes much due to the gzip compression algorithm used in ROOT to save equal-value variables only once.
In principle, electric field values can be calculated from the potential using Eq. 4 when needed. However, given their frequent usage, their values on each point are calculated and saved in GeFiCa after the SOR calculation for the potential.
For a three dimensional grid with a variable step size, 14 std::vector s with a double precision are needed in total to save 3 coordinates, 3 × 2 distances to previous and next points, 1 potential, 1 total electric field and its 3 components. They are public member variables in the class Points inherited by all grid classes shown in Fig. 24, including those representing lower dimensional grids, where variables for higher dimensions are not used at all. Since the C++ vector does not allocate memory if it has no element, there is no penalty in storage size in this solution. An alternative is to create Points1D, Points2D, Point3D, and consequently, Grid1D, Grid2D, Grid3D for various dimensions. This complicates the overall class structure unnecessarily, hence is not used in GeFiCa.
A few more vector s are added in the Grid class to record space charge densities in individual grid points, as well as flags to tell whether a point is in or out of a crystal, and whether it is in or out of the depleted region.
As described in Section 5.5, an electric field line can be saved in a series of points with variable distances between them. That is why the class FieldLine inherits the data structure from Points, as shown in Fig. 24.
The Grid and the Detector classes are both daughters of the TNamed class in ROOT, which inherits the capability to stream its data members for I/O from the TObject class in ROOT. Consequently, all concrete grid and detector classes can be directly saved into a standard ROOT file using one line of C++ as shown in the following code snippet: The first line creates a TTree object t out of the saved field values in the rhoz object. The second line draws the potential, v, on the first (y-axis) and second (x-axis) coordinates, c1 and c2, as a colored contour along zaxis (the "colz" option), as shown in Fig. 25.
The price to pay for all these convenience is that the objects saved in the ROOT file can only be loaded without warning message when the compiled GeFiCa library can be found and automatically loaded by ROOT. The way to realize this is detailed in Section. 6.6.

Detector Configurations
Two pieces of information are needed for electric field calculation: first, boundary conditions, and second, the space charge distribution. Fig. 25 The color contour of the potential field of an ICPC detector drawn with TTree::Draw("c2:c1:v","","colz").
Boundary conditions can be set through the detector geometry and voltages on electrodes. Take the previously defined PointContact detector as an example, its basic dimensions can be set as In case of a segmented detector, the bias voltage array can have more than two elements. The index of an element can be kept the same as the corresponding segment identification number.
As described in detail in Sec. 2, it is reasonable to use a first-order polynomial to approximate the space charge distribution in a HPGe crystal. With this simplification, we just need to specify the impurities at the top and the bottom of a crystal given by the manufacturer. For example, d e t e c t o r . BottomImpurity=3e9 /cm3 ; d e t e c t o r . TopImpurity=7e9 /cm3 ; The impurity level at a specific axial position is interpolated in GeFiCa based on these two numbers. In case of a small crystal, the impurity can be regarded as a constant. Its average impurity can be set as d e t e c t o r . S e t A v e r a g e I m p u r i t y ( 3 e9 /cm3 ) ; Fig. 26 Cross section along z-axis of an inverted coaxial pointcontact HPGe detector, and parameters describing its dimensions. To shorten the names, width is represented as a single capital case "W", height "H", and radius "R".

Units and Constants
We have seen in previous code snippets that an input parameter in GeFiCa is a product of a number and a unit. Common units and constants for field calculation, together with their conversion rules, are defined in GeFiCa/src/Units.h. The following is a snippet of the file: The advantage of this unit system is three-fold. First, the code is self-explainable, there is no ambiguity in the unit of an input value. Second, the user has freedom to choose units, such as "mm" instead of "cm", or "kV" instead of "volt". Otherwise, he or she has to use the set of units used for internal calculation. Third, since the unit conversion rules are pre-defined, there is no need to worry about them when using input parameters for internal calculations. The programmer can focus on the logic instead of unit conversion. This way of handling units is adopted from Geant4 [22][23][24]. Most of the units and constants have been defined in Geant4 already. However, since only a small subset of the units are useful for field calculation, they are re-defined in GeFiCa to avoid unnecessary dependence on Geant4.

Compilation and Installation
GeFiCa relies on ROOT to realize C++ and Python scripting, efficient I/O and plotting. It has to be compiled against ROOT libraries. This is achieved through a simple Makefile that uses the root−config executable available from any successfully installed ROOT package to get the location of ROOT libraries and necessary compilation flags. The compilation process is as simple as cd / path / t o /GeFiCa/ s r c && make After a successful compilation a shared C++ library, libGeFiCa.so, can be found in /path/to/GeFiCa/src. Once its location is added to the LD_LIBRARY_PATH environment variable (or DYLD_LIBRARY_PATH in MacOS), the library can be automatically loaded only when needed as any other ROOT libraries in ROOT interactive sessions or scripts thanks to the rootmap and pcm files [44] generated by the make process in the same directory as the library.

Supported OS
Since ROOT is available in the three common operating systems, Linux, Windows and MacOS, in principle, GeFiCa should be able to be compiled in all of them as well. However, since GeFiCa relies on a simple Makefile to compile, it cannot be directly compiled through the native Windows compilation system. Instead, it can be compiled in a Windows Subsystem for Linux (WSL). To date, GeFiCa has been compiled successfully in CentOS 6 and 7, MacOS 10.12, 10.13, 10.14, and Ubuntu 18.04 as a WSL.

Code Accessibility
The codes of GeFiCa are hosted online at GitHub [45]: https://github.com/jintonic/gefica. They can be downloaded directly from the web page or through git. GeFiCa is release under the MIT License [46]. It can be freely used without any warranty as long as the license is distributed along with it.

Code Documentation
A git branch gh-pages is used to host the homepage code for GeFiCa. The homepage is available under a cus-tomized domain name: http://physino.xyz/gefica. It lists three main resources about GeFiCa that one can get help from: the GeFiCa repository page hosted on GitHub, the code documentation hosted on https:// codedocs.xyz, and the user manual hosted on https: //readthedocs.org.
There is a README.md file in each directory in GeFiCa to explain the contents of the directory written in GitHub Flavored Markdown format [47]. They are rendered to web pages automatically in GitHub. A user can quickly get help with or without the source code.
Explanations of GeFiCa classes and variables are embedded in the source code as C++ comments using the Doxygen [48] convention. They can be rendered by Doxygen into nicely formatted documentations locally or on https://codedocs.xyz. The online version is updated automatically once new codes are pushed to the GitHub repository.
The user manual is written in restructured text format and can be rendered to web pages locally or on https://readthedocs.org. The online version is updated automatically once new documentation is pushed to the GitHub repository.
In addition to these, example codes are shipped with GeFiCa as ROOT macros as described in detail in the next section to demonstrate the usage of individual GeFiCa classes.

Macros and Scripts
A modern C++ interpreter, cling [49], has been created and adopted as the back-end of the interactive session of ROOT [33] since the version 6 of it. A user can run C++ snippets, sometimes called ROOT macros or scripts, interactively in cling without writing and compiling the "main" function. With immediate feedback after the execution of each line of a script, a user can learn and experiment a new C++ class, a function, or simply a syntax easily. To fulfill its educational purpose, GeFiCa is compiled as a ROOT library. All snippets in previous sections demonstrating the configuration of a detector or the operation of a grid can be run as they are in cling.
ROOT also provides a Python extension module, PyROOT, that allows the user to interact with any ROOT class from the Python interpreter. For users who prefer the Python interpreter to cling, they can call GeFiCa classes with Python syntax directly in the standard Python interpreter.
It is worth noting that cling comes with a Jupyter [50] kernel, which makes it possible to run GeFiCa scripts in a Jupyter notebook with either C++ or python syntax.
All concrete grid and detector classes in GeFiCa inherit the capability to inspect themselves from the TObject class in ROOT. Some standard functions in TObject, such as Dump(), can be used to check the default or user-specified configurations of a grid or detector object, as shown in Listing 1. The first column of the output are the member variables of the GeFiCa::X class. The second are their current values. The last are explanations of those variable. These explanations are written as C++ comments after the member variables. They can be parsed by both Doxygen and ROOT to generate code documentation in various formats and contexts.
Macros are organized in sub-folders in GeFiCa/examples/ to demonstrate the usage of GeFiCa classes. The planar/, trueCoaxial/, hemispherical/, pointContact/, and segmented/ folders are used to show how to configure specific types of HPGe detectors and then calculate the fields in them. The analytic/ and the fenics/ folders contains macros that are independent of the GeFiCa libraries. The macros in the former demonstrate how to calculate and visualize the field distribution in simple HPGe detectors using ROOT. The latter shows Python codes to calculate and visualize the field distribution in a simplified point-contact geometry using FEniCS [19]. All field distributions shown in this work are generated using these macros. A user can learn the topics by at first running these macros to reproduce plots in this work, and then modifying them to meet his/her own needs.

Code Verification
A common way to verify the saneness of a complex theory in physics is to consider extreme conditions, under which the theory can be simplified and compared to predictions based on common sense. Take the field in a point-contact detector as an example, there are two extreme cases where the field in certain part of the detector can be regarded as the same as that in a planar or a true-coaxial detector. This makes it possible to compare the numeric calculation of a point-contact detector field directly with analytic solutions.

Comparison with Analytic Solutions
In the first extreme case we consider a point-contact detector that takes a pancake-like shape, that is, its thickness is much smaller than its diameter. Furthermore, its "point-contact" covers almost the entire bottom end surface. The electric potential in such a detector is shown in the bottom plot in Fig. 27. At the radial center of the detector, the field is essentially the same as that in a planar detector that has the same thickness and impurity concentration. In the top plot in Fig. 27 the analytic solution of such a planar detector is overlaid on top of the numerical result of the pancake-like "point-contact" detector along the axial positions at radius = 0. Fig. 27 Top: Comparison of the electric potential calculated numerically in a pancake-like "point-contact" detector with the analytic solution of a planar detector that has the same thickness (or height) and impurity concentration. Bottom: The electric potential distribution calculated numerically in the pancake-like detector, the "point-contact" of which is artificially enlarged to cover almost the entire bottom end surface.
In the second case let us consider a point-contact detector that looks like a thin stick, that is, its diameter is much smaller than its height. Furthermore, let's make its "point-contact" as deep into the crystal as possible. The potential distribution in such a detector calculated numerically is shown in the right plot in Fig. 28. Far away from the top end of the detector, the field is essentially the same as that in a true-coaxial detector that has the same radius and impurity concentration. In the left plot in Fig. 28 the analytic solution of such a true-coaxial detector is overlaid on top of the numerical result of the thin-stick-like "point-contact" detector along the radius at an axial position 5 mm above the bottom surface.

Comparison with Fieldgen
Even though the perfect matches between the analytic solutions and the numerical results in both cases are convincing evidences of the correctness of the numerical calculation implemented in GeFiCa, it is worth noting that a constant impurity concentration throughout the entire crystal is assumed to make the analytic solutions possible. In case of an arbitrary impurity dis-  Left: Comparison of the electric potential calculated numerically in a thin-stick-like "point-contact" detector with the analytic solution of a true-coaxial detector that has the same radius and impurity concentration. Right: The electric potential distribution calculated numerically in the thin-stick-like detector, the "point-contact" of which is artificially prolonged along almost the entire height of the crystal. tribution, no simple analytic solution is available, the numerical calculation in GeFiCa is compared to that of fieldgen [1,14], a thoroughly examined and widely accepted package in the field, given identical point-contact detector configurations.
The biggest difference between GeFiCa and fieldgen in the aspect of numerical calculation is probably the setup of grid points. In case of fieldgen, the grid points along the radial direction, r, of a detector start from r = 0 and end at r = the radius of the detector. In case of GeFiCa, the grid points are in the range of [−radius, +radius] and there is no grid point at r = 0 to avoid setting artificial boundary conditions at r = 0. Due to this difference, there is no one-to-one correspondence between a grid point in GeFiCa and a grid point in fieldgen. In order to make a point-to-point comparison, linear interpolation is used to get the total electric field strength at a fieldgen point from two nearby GeFiCa points, the interpolated value is then compared to the fieldgen value at the same point. Their relative difference in percentage is shown as colored contour in Fig. 29.
The largest difference is about 8.5% at the top right corner of the point-contact. This point is removed from Fig. 29 so that subtle differences between fieldgen and GeFiCa are more visible in the figure. The second largest difference is about 2.5% at an adjacent point, shown as the red spot in Fig. 29. The difference quickly falls below 0.1% only a few points away from the corner, which translates to about one mini-meter in length given the 0.1 mm distance between grid points. Such difference is most probably due to different treatments in fieldgen and GeFiCa on grid points near boundaries.
Fortunately, the difference is of little importance in practice since there is no such sharp corner inside any detector in reality. Predictions of GeFiCa and fieldgen in the bulk of the detector are essentially identical.

Relaxation Factor
As described in Sec. 5.2, the number of iterations needed for a successive relaxation process to converge can be reduced by introducing a relaxation factor in between [1,2). Fig. 30 shows the number of iterations for a succesive over-relaxation (SOR) process to converge as a function of the relaxation factor. Each data point in the figure represents the result from a numerical calculation of the field in an ideal planar detector. Data points that  A common trend shared by all the lines is that there is a point where the number of iterations is minimized. That is where GeFiCa reaches its best performance. As the number of grid points increase from 101 to 601, the relaxation factor corresponding to the minimal iteration numbers increases from around 1.94 to 1.99. The default value of the relaxation factor is set to 1.95 in GeFiCa. A user can change it using the following line of code if desired. g r i d . R e l a x a t i o n F a c t o r = 1 . 9 9 ; The gain in performance by selecting an appropriate relaxation factor becomes more prominent when the number of grid points becomes larger. Take the up most curve in Fig. 30 as an example, which corresponds to calculations done with the finest grid, when the relaxation factor changes by only 0.06 from 1.93 to 1.99, the number of iterations reduces from about 10,000 to less than 2,000. While the lowest curve in Fig. 30 is almost flat around 1.94, that is, the relaxation factor cannot help much to gain speed for calculations with a very coarse grid. This is not a problem since those calculations are fast already.
After every 100 iterations, GeFiCa prints the overall difference of potentials at all grid points between current and previous iterations. When the difference is smaller than a target precision (1 × 10 −7 V by default), the SOR is regarded as converged, the calculation stops there, and the CPU time used for the calculatioin is printed out on screen as shown in the terminal output below: This terminal output is associated with the calculation used to generate Fig. 29. The overall number of grid points is 349,140. The CPU time used for the calculation is about 23 second in a Linux server with an Intel Xeon Gold 5118 CPU at 1 GHz. The relaxation factor chosen for this calculation is 1.994.

Output File Size
The output of fieldgen used to generate Fig. 29 is saved as a simple ASCII file that is 8.1 Mega bytes in size. The detector configuration is saved as a short header of the file. The rest of the file are six columns of values of the grid point positions (radial and axial), the voltage, the overall electric field strength, and its radial and axial components.
As described in Sec. 6.3, the detector and grid objects in GeFiCa can be directly saved in a standard ROOT file. Its contents can be printed and visualized in a ROOT interactive session as demonstrated in the code snippets in Sec. 6.3 and 6.10. In addition to the information saved in a fieldgen output, a GeFiCa output also contains the intervals between grid points, flags indicating whether a point is depleted or not, etc. It also contains about twice more grid points than fieldgen. In total, the amount of information saved in GeFiCa is about 4 times more than that saved in the fieldgen output. The size of the GeFiCa output ROOT file used to generate Fig. 29 is 9.2 Mega bytes, only slightly larger than that of the fieldgen ouput file, thanks to the gzip algrithm used to compress a ROOT file mentioned in Sec. 6.3.

Extendability and Limitation
Let's take a realistic planar detector configuration shown in Fig. 31 as an example to demonstrate the procedure of extending GeFiCa for a new type of detector. The top and bottom surfaces of the detector are covered with a thin layer of aluminium to form the electric contacts. All the side surfaces are covered with a thin layer of amorphous germanium for passivation purpose. The two side wings can be used for handling the detector without touching its sensitive surfaces [51]. Since they are thin compared to the overall thickness of the detector, the electric field distribtuion inside the detector can hence be approximated by that in an ideal 1D planar detector. However, if our intention is to study the influence of the thickness of the wings on the electric field, we need at least a 2D grid in Cartesian coordinates to perform the numerical calculation, which can be achieved with the following steps.
At first, a class called XY that represents the dimension and coordinates needs to be created. It inherits all member variables in its base class Grid that define the grid. Since the numerical expression of Poisson's Equation (Eq. 5) depends on dimensions and coordinates used for the calculation, a protected virtual function, void OverRelaxAt(size_t idx), in Grid needs to be overwritten in XY, which takes care of the updating of the field value at each grid point indexed by idx.
Secondly, a class called TopHat that describes the geometry of the detector needs to be created. It inherits the member variables that hold voltages values of all electrodes from its base class Detector. It also inherits the impurity distribution from the class, crystal. A public member function void Draw() in Detector needs to be overwritten in TopHat to visualize the geometry setup.
At last the public virtual function in Grid called void SetupWith(Detector&) needs to be overwritten in XY, which takes the boundary conditions and impurity distribution from TopHat to construct and initialize the grid for the calculation.
For completeness, a folder called TopHat is recommended to be created under GeFiCa/examples, which contains ROOT or Python scripts demonstrating the usage of XY and TopHat.
Given its extendability, there is no limitation on GeFiCa from the functionality point of view. From the education point of view, however, there is currently no function in GeFiCa demonstrating the adaptive grid configuration that automatically updates distances between grid points over iterations based on the strength of local electric field. Note that there is no fundamental limitation from GeFiCa inhibiting doing so, since there are separated member variables in the grid class to hold distances from a grid point to its neighbors in all directions. Practically, GeFiCa is already fast and precise enough with fixed step length for common HPGe configurations. This function can be added in if necessary.

Summary
The new educational program, GeFiCa, has been created to demonstrate analytic and numeric methods to calculate static electric fields and potentials in HPGe detectors. It is freely available from http://physino. xyz/gefica and can be installed in three major operating systems, Linux, MacOS and Windows, as a CERN ROOT [33] library extension. Powered by ROOT, GeFiCa allows its users to explore in detail the calculation procedure by executing C++ or Python code snippets in ROOT interactive sessions or Jupyter notebooks without compilation. Example code snippets are shipped together with the library to demonstrate calculations for common detector configurations, and to visualize the resulting field distributions in graphs or color contours. In addition to field calculations, GeFiCa offers functionalities to calculate the HPGe detector depletion voltage, undepleted region, capacitance, etc., that are not available from general-purpose field calculation programs, such as Maxwell3D and FEniCS. Compared to open projects that are also specialized in HPGe field calculation, such as fieldgen and SSD, etc., GeFiCa offers a ROOT-based C++ solution that is equally accurate and efficient, and shipped with a large amount of documentations and examples that are not readily available in others.
This article was written to provide an entry level review of methods and tools available at the moment, with the hope that its readers feel comfortable to make an educated choice of simulation tools best suited for the task at their hands.
where x, y, z are coordinates, dx + , dy + , dz + are distances to the next grid points, dx − , dy − , dz − distances to the previous. In a 1D cylindrical coordinate, the potential at a grid point after the i-th successive relaxation iteration, V i+1 , can be expressed as where r is the coordinate. dr + is the distance to the next grid point, dr − the distance to the previous. In 3D cylindrical coordinate, the potential at a grid point after the i-th successive relaxation iteration, V i+1 , can be expressed as where r, θ, z are coordinates, dr + , dθ + , dz + , are step lengths to the next grid points, dr − , dθ − , dz − step lengths to the previous. In a 1D spherical coordinate, the potential at a grid point after the i-th successive relaxation iteration, V i+1 , can be expressed as where r is the coordinate, dr + is the distance to the next grid point, dr − the distance to the previous. In 3D Spherical Coordinate, the potential at a grid point after the i-th successive relaxation iteration, V i+1 , can be expressed as where r, θ, φ are coordinates, dr + , dθ + , dφ + , are the step lengths to the next grid points, dr − , dθ − , dφ − to the previous.