Merge branch 'staging' into 'master'

* Andreas Matthias Nürnberg, KIT, [nurnberg](https://gitlab.cern.ch/nurnberg)

* Sebastian Pape, TU Dortmund University, [spape](https://gitlab.cern.ch/spape)

* Marko Petric, CERN, [mpetric](https://gitlab.cern.ch/mpetric)

* Florian Michael Pitters, HEPHY, [fpipper](https://gitlab.cern.ch/fpipper)

* Radek Privara, Palacky University Olomouc, [rprivara](https://gitlab.cern.ch/rprivara)

* Nashad Rahman, The Ohio State University, [nashadroid](https://github.com/nashadroid)

* Sabita Rao, GSDocs2020 Student, [srao](https://gitlab.cern.ch/srao)

url = {https://www.sciencedirect.com/science/article/pii/088328899190002I},

author = {S. Croft and D.S. Bond}

}

@article{massey,

  author  = {D. J. {Massey} and J. P. R. {David} and G. J. {Rees}},

  journal = {IEEE Transactions on Electron Devices},

  title   = {Temperature Dependence of Impact Ionization in Submicrometer Silicon Devices},

  year    = {2006},

  volume  = {53},

  number  = {9},

  pages   = {2328-2334},

  doi     = {10.1109/TED.2006.881010},

  ISSN    = {1557-9646}

}

@article{overstraeten,

  title = {Measurement of the ionization rates in diffused silicon p-n junctions},

  journal = {Solid-State Electronics},

  volume = {13},

  number = {5},

  pages = {583-608},

  year = {1970},

  issn = {0038-1101},

  doi = {10.1016/0038-1101(70)90139-5},

  url = {https://www.sciencedirect.com/science/article/pii/0038110170901395},

  author = {R. {Van Overstraeten} and H. {De Man}},

}

@article{okuto,

  title = {Threshold energy effect on avalanche breakdown voltage in semiconductor junctions},

  journal = {Solid-State Electronics},

  volume = {18},

  number = {2},

  pages = {161-168},

  year = {1975},

  issn = {0038-1101},

  doi = {10.1016/0038-1101(75)90099-4},

  url = {https://www.sciencedirect.com/science/article/pii/0038110175900994},

  author = {Y. Okuto and C.R. Crowell},

}

@inproceedings{bologna,

  author={Valdinoci, M. and Ventura, D. and Vecchi, M.C. and Rudan, M. and Baccarani, G. and Illien, F. and Stricker, A. and Zullino, L.},

  booktitle={1999 International Conference on Simulation of Semiconductor Processes and Devices. SISPAD'99 (IEEE Cat. No.99TH8387)},

  title={Impact-ionization in silicon at large operating temperature},

  year={1999},

  volume={},

  number={},

  pages={27-30},

  doi={10.1109/SISPAD.1999.799251}

 }

---

# SPDX-FileCopyrightText: 2022 CERN and the Allpix Squared authors

# SPDX-License-Identifier: CC-BY-4.0

title: "Impact Ionization"

weight: 5

---

Allpix Squared implements charge multiplication via impact ionization models. These models are only used by propagation

modules which perform a step-by-step simulation of the charge carrier motion.

The gain $`g`$ is calculated for all models as exponential of the model-dependent impact ionization coefficient $`\alpha`$ and

the length of the step $`l`$ performed in the respective electric field. If the electric field strength stays below a

configurable threshold $`E_{\text{thr}}`$, unity gain is assumed:

```math

g (E, T) = \left\{

\begin{array}{ll}

    e^{l \cdot \alpha(E, T)} & E > E_{\text{thr}}\\

    1.0 & E < E_{\text{thr}}

\end{array}

\right.

```

The following impact ionization models are available:

## Massey Model

The Massey model \[[@massey]\] describes impact ionization as a function of the electric field $`E`$.

The ionization coefficients are parametrized as

```math

\alpha (E, T) = A e^{-\frac{B(T)}{E}},

```

where $`A`$ and $`B(T)`$ are phenomenological parameters, defined for electrons and holes respectively.

While $`A`$ is assumed to be temperature-independent, parameter $`B`$ exhibits a temperature dependence and is defined as

```math

B(T) = C + D \cdot T.

```

The parameter values implemented in Allpix Squared are taken from Section 3 of \[[@massey]\] as:

```math

\begin{aligned}

    A_{e} &= 4.43\times 10^{5} \,\text{/cm}\\

    C_{e} &= 9.66\times 10^{5} \,\text{V/cm}\\

    D_{e} &= 4.99\times 10^{2} \,\text{V/cm/K}\\

\\

    A_{h} &= 1.13\times 10^{6} \,\text{/cm}\\

    C_{h} &= 1.71\times 10^{6} \,\text{V/cm}\\

    D_{h} &= 1.09\times 10^{3} \,\text{V/cm/K}\\

\end{aligned}

```

for electrons and holes, respectively.

This model can be selected in the configuration file via the parameter `multiplication_model = "massey"`.

## Van Overstraeten-De Man Model

The Van Overstraeten-De Man model \[[@overstraeten]\] describes impact ionization using Chynoweth's law, given by

```math

\alpha (E, T) = \gamma (T) \cdot a_{\infty} \cdot e^{-\frac{\gamma(T) \cdot b}{E}},

```

For holes, two sets of impact ionization parameters $`p = \left\{ a_{\infty}, b \right\}`$ are used depending on the electric field:

```math

p = \left\{

\begin{array}{ll}

    p_{\text{low}} & E < E_{0}\\

    p_{\text{high}} & E > E_{0}

\end{array}

\right.

```

Temperature scaling of the ionization coefficient is performed via the $`\gamma(T)`$ parameter following the Synposys

Sentaurus TCAD user manual as:

```math

\gamma (T) = \tanh \left(\frac{0.063\times 10^{6} \,\text{eV}}{2 \times 8.6173\times 10^{-5} \,\text{eV/K} \cdot T_0} \right) \cdot \tanh \left(\frac{0.063\times 10^{6} \,\text{eV}}{2 \times 8.6173\times 10^{-5} \,\text{eV/K} \cdot T} \right)^{-1}

```

The value of the reference temperature $`T_0`$ is not entirely clear as it is never stated explicitly, a value of

$`T_0 = 300 \,\text{K}`$ is assumed. The other parameter values implemented in Allpix Squared are taken from the abstract

of \[[@overstraeten]\] as:

```math

\begin{aligned}

    E_0 &= 4.0\times 10^{5} \,\text{V/cm}\\

    a_{\infty, e} &= 7.03\times 10^{5} \,\text{/cm}\\

    b_{e} &= 1.231\times 10^{6} \,\text{V/cm}\\

    \\

    a_{\infty, h, \text{low}} &= 1.582\times 10^{6} \,\text{/cm}\\

    a_{\infty, h, \text{high}} &= 6.71\times 10^{5} \,\text{/cm}\\

    b_{h, \text{low}} &= 2.036\times 10^{6} \,\text{V/cm}\\

    b_{h, \text{high}} &= 1.693\times 10^{6} \,\text{V/cm}

\end{aligned}

```

This model can be selected in the configuration file via the parameter `multiplication_model = "overstraeten"`.

## Okuto-Crowell Model

The Okuto-Crowell model \[[@okuto]\] defines the impact ionization coefficient similarly to the above models but in addition

features a linear dependence on the electric field strength $`E`$. The coefficient is given by:

```math

\alpha (E, T) = a(T) \cdot E \cdot e^{-\left(\frac{b(T)}{E}\right)^2}.

```

The two parameters $`a, b`$ are temperature dependent and scale with respect to the reference temperature

$`T_0 = 300 \,\text{K}`$ as:

```math

\begin{aligned}

    a(T) &= a_{300} \left[ 1 + c\left(T - T_0\right) \right]\\

    b(T) &= a_{300} \left[ 1 + d\left(T - T_0\right) \right]

\end{aligned}

```

The parameter values implemented in Allpix Squared are taken from Table 1 of \[[@okuto]\], using the values for silicon, as:

```math

\begin{aligned}

    a_{300, e} &= 0.426 \,\text{/V}\\

    c_{e} &= 3.05\times 10^{-4}\\

    b_{300, e} &= 4.81\times 10^{5} \,\text{V/cm}\\

    d_{e} &= 6.86\times 10^{-4}\\

\\

    a_{300, h} &= 0.243 \,\text{/cm}\\

    c_{h} &= 5.35\times 10^{-4}\\

    b_{300, h} &= 6.53\times 10^{5} \,\text{V/cm}\\

    d_{h} &= 5.67\times 10^{-4}\\

\end{aligned}

```

This model can be selected in the configuration file via the parameter `multiplication_model = "okuto"`.

## Bologna Model

The Bologna model \[[@bologna]\] describes impact ionization for experimental data in an electric field range from

$`130 \,\text{kV/cm}`$ to $`230 \,\text{kV/cm}`$ and temperatures up to $`400 \,^{\circ}\text{C}`$. The impact ionization

coefficient takes a different form than the previous models and is given by

```math

\alpha (E, T) = \frac{E}{a(T) + b(T) e^{d(T) / \left(E + c(T) \right)}},

```

for both electrons and holes.

The temperature-dependent parameters $`a(T), b(T), c(T)`$ and $`d(T)`$ are defined as:

```math

\begin{aligned}

    a(T) &= a_{0} + a_1 T^{a_2}\\

    b(T) &= b_{0} e^{b_1 T}\\

    c(T) &= c_{0} + c_1 T^{c_2} + c_3 T^{2}\\

    d(T) &= d_{0} + d_1 T + d_2 T^{2}\\

\end{aligned}

```

The parameter values implemented in Allpix Squared are taken from Table 1 of \[[@bologna]\] as:

```math

\begin{aligned}

    a_{0, e} &= 4.3383 \,\text{V}\\

    a_{1, e} &= -2.42\times 10^{-12} \,\text{V}\\

    a_{2, e} &= 4.1233\\

    b_{0, e} &= 0.235 \,\text{V}\\

    b_{1, e} &= 0\\

    c_{0, e} &= 1.6831\times 10^{4} \,\text{V/cm}\\

    c_{1, e} &= 4.3796 \,\text{V/cm}\\

    c_{2, e} &= 1\\

    c_{3, e} &= 0.13005 \,\text{V/cm}\\

    d_{0, e} &= 1.2337\times 10^{6} \,\text{V/cm}\\

    d_{1, e} &= 1.2039\times 10^{3} \,\text{V/cm}\\

    d_{2, e} &= 0.56703 \,\text{V/cm}\\

\\

    a_{0, h} &= 2.376 \,\text{V}\\

    a_{1, h} &= 1.033\times 10^{-2} \,\text{V}\\

    a_{2, h} &= 1\\

    b_{0, h} &= 0.17714 \,\text{V}\\

    b_{1, h} &= -2.178\times 10^{-3} \,\text{/K}\\

    c_{1, h} &= 0\\

    c_{1, h} &= 9.47\times 10^{-3} \,\text{V/cm}\\

    c_{2, h} &= 2.4924\\

    c_{3, h} &= 0\\

    d_{0, h} &= 1.4043\times 10^{6} \,\text{V/cm}\\

    d_{1, h} &= 2.9744\times 10^{3} \,\text{V/cm}\\

    d_{2, h} &= 1.4829 \,\text{V/cm}\\

\end{aligned}

```

This model can be selected in the configuration file via the parameter `multiplication_model = "bologna"`.

[@massey]: https://doi.org/10.1109/TED.2006.881010

[@overstraeten]: https://doi.org/10.1016/0038-1101(70)90139-5

[@okuto]: https://doi.org/10.1016/0038-1101(75)90099-4

[@bologna]: https://doi.org/10.1109/SISPAD.1999.799251

    config_.setDefault<bool>("ignore_magnetic_field", false);

    // Set defaults for charge carrier multiplication

    config_.setDefault<std::string>("multiplication_model", "none");

    config_.setDefault<double>("multiplication_threshold", 1e-2);

    // Copy some variables from configuration to avoid lookups:

    temperature_ = config_.get<double>("temperature");

    timestep_min_ = config_.get<double>("timestep_min");

                                                     static_cast<int>(Units::convert(integration_time_, "ns") * 5),

                                                     0,

                                                     static_cast<double>(Units::convert(integration_time_, "ns")));

        gain_histo_ = CreateHistogram<TH1D>(

            "gain_histo", "Gain per charge carrier group after propagation;gain;number of groups transported", 500, 1, 25);

    }

    // Prepare mobility model

    // Prepare recombination model

    recombination_ = Recombination(config_, detector_->hasDopingProfile());

    // Impact ionization model

    multiplication_ = ImpactIonization(config_);

    // Check multiplication and step size larger than a picosecond:

    if(!multiplication_.is<NoImpactIonization>() && timestep_max_ > 0.001) {

        LOG(WARNING) << "Charge multiplication enabled with maximum timestep larger than 1ps" << std::endl

                     << "This might lead to unphysical gain values.";

    }

    // Prepare trapping model

    trapping_ = Trapping(config_);

}

            }

            // Propagate a single charge deposit

            auto [final_position, time, state] = propagate(

            auto [final_position, time, gain, state] = propagate(

                initial_position, deposit.getType(), deposit.getLocalTime(), event->getRandomEngine(), output_plot_points);

            if(state == CarrierState::RECOMBINED) {

            }

            LOG(DEBUG) << " Propagated " << charge_per_step << " to " << Units::display(final_position, {"mm", "um"})

                       << " in " << Units::display(time, "ns") << " time, final state: " << allpix::to_string(state);

                       << " in " << Units::display(time, "ns") << " time, gain " << gain

                       << ", final state: " << allpix::to_string(state);

            // Create a new propagated charge and add it to the list

            auto global_position = detector_->getGlobalPosition(final_position);

            PropagatedCharge propagated_charge(final_position,

                                               global_position,

                                               deposit.getType(),

                                               charge_per_step,

                                               static_cast<unsigned int>(std::round(charge_per_step * gain)),

                                               deposit.getLocalTime() + time,

                                               deposit.getGlobalTime() + time,

                                               state,

 * velocity at every point with help of the electric field map of the detector. An Runge-Kutta integration is applied in

 * multiple steps, adding a random diffusion to the propagating charge every step.

 */

std::tuple<ROOT::Math::XYZPoint, double, CarrierState>

std::tuple<ROOT::Math::XYZPoint, double, double, CarrierState>

GenericPropagationModule::propagate(const ROOT::Math::XYZPoint& pos,

                                    const CarrierType& type,

                                    const double initial_time,

    // Create a runge kutta solver using the electric field as step function

    Eigen::Vector3d position(pos.x(), pos.y(), pos.z());

    // Initialize gain

    double gain = 1.;

    // Define a function to compute the diffusion

    auto carrier_diffusion = [&](double efield_mag, double doping_concentration, double timestep) -> Eigen::Vector3d {

        double diffusion_constant = boltzmann_kT_ * mobility_(type, efield_mag, doping_concentration);

    // Continue propagation until the deposit is outside the sensor

    Eigen::Vector3d last_position = position;

    ROOT::Math::XYZVector efield{}, last_efield{};

    double last_time = 0;

    size_t next_idx = 0;

    auto state = CarrierState::MOTION;

        // Save previous position and time

        last_position = position;

        last_time = runge_kutta.getTime();

        last_efield = efield;

        // Execute a Runge Kutta step

        auto step = runge_kutta.step();

        // Get the current result and timestep

        auto timestep = runge_kutta.getTimeStep();

        position = runge_kutta.getValue();

        LOG(TRACE) << "Step from " << Units::display(static_cast<ROOT::Math::XYZPoint>(last_position), {"um"}) << " to "

                   << Units::display(static_cast<ROOT::Math::XYZPoint>(position), {"um"});

        // Get electric field at current position and fall back to empty field if it does not exist

        auto efield = detector_->getElectricField(static_cast<ROOT::Math::XYZPoint>(position));

        efield = detector_->getElectricField(static_cast<ROOT::Math::XYZPoint>(position));

        auto doping = detector_->getDopingConcentration(static_cast<ROOT::Math::XYZPoint>(position));

        // Apply diffusion step

                   << " to " << Units::display(static_cast<ROOT::Math::XYZPoint>(position), {"um", "mm"}) << " at "

                   << Units::display(initial_time + runge_kutta.getTime(), {"ps", "ns", "us"})

                   << ", state: " << allpix::to_string(state);

        // Adapt step size to match target precision

        double uncertainty = step.error.norm();

        // Apply multiplication step, fully deterministic from local efield and step length; Interpolate efield values

        gain *= multiplication_(type, (std::sqrt(efield.Mag2()) + std::sqrt(last_efield.Mag2())) / 2., step.value.norm());

        if(gain > 20.) {

            LOG(WARNING) << "Detected gain of " << gain << ", local electric field of "

                         << Units::display(std::sqrt(efield.Mag2()), "kV/cm") << ", diode seems to be in breakdown";

        } else if(gain > 1.) {

            LOG(DEBUG) << "Calculated gain of " << gain << " for step of " << Units::display(step.value.norm(), {"um", "nm"})

                       << " from field of " << Units::display(std::sqrt(last_efield.Mag2()), "kV/cm") << " to "

                       << Units::display(std::sqrt(efield.Mag2()), "kV/cm");

        }

        // Update step length histogram

        if(output_plots_) {

            uncertainty_histo_->Fill(static_cast<double>(Units::convert(step.error.norm(), "nm")));

        }

        // Adapt step size to match target precision

        double uncertainty = step.error.norm();

        // Lower timestep when reaching the sensor edge

        if(std::fabs(model_->getSensorSize().z() / 2.0 - position.z()) < 2 * step.value.z()) {

            timestep *= 0.75;

        }

    }

    if(output_plots_) {

        gain_histo_->Fill(gain);

    }

    if(state == CarrierState::RECOMBINED) {

        LOG(DEBUG) << "Charge carrier recombined after " << Units::display(last_time, {"ns"});

    } else if(state == CarrierState::TRAPPED) {

    }

    // Return the final position of the propagated charge

    return std::make_tuple(static_cast<ROOT::Math::XYZPoint>(position), initial_time + time, state);

    return std::make_tuple(static_cast<ROOT::Math::XYZPoint>(position), initial_time + time, gain, state);

}

void GenericPropagationModule::finalize() {

        trapped_histo_->Write();

        recombination_time_histo_->Write();

        trapping_time_histo_->Write();

        gain_histo_->Write();

    }

    long double average_time = static_cast<long double>(total_time_picoseconds_) / 1e3 /

#include "objects/DepositedCharge.hpp"

#include "objects/PropagatedCharge.hpp"

#include "physics/ImpactIonization.hpp"

#include "physics/Mobility.hpp"

#include "physics/Recombination.hpp"

#include "physics/Trapping.hpp"

         * @param initial_time Initial time passed before propagation starts in local time coordinates

         * @param random_generator Reference to the random number engine to be used

         * @param output_plot_points Reference to vector to hold points for line graph output plots

         * @return Tuple with the point where the deposit ended after propagation, the time the propagation took and the

         * final state of the charge carrier at the end of processing

         * @return Tuple with the point where the deposit ended after propagation, the time the propagation took, the

         * cumulative gain and the final state of the charge carrier at the end of processing

         */

        std::tuple<ROOT::Math::XYZPoint, double, CarrierState> propagate(const ROOT::Math::XYZPoint& pos,

        std::tuple<ROOT::Math::XYZPoint, double, double, CarrierState> propagate(const ROOT::Math::XYZPoint& pos,

                                                                                 const CarrierType& type,

                                                                                 const double initial_time,

                                                                                 RandomNumberGenerator& random_generator,

        // Models for electron and hole mobility and lifetime

        Mobility mobility_;

        Recombination recombination_;

        ImpactIonization multiplication_;

        Trapping trapping_;

        // Precalculated value for Boltzmann constant:

        Histogram<TH1D> trapped_histo_;

        Histogram<TH1D> recombination_time_histo_;

        Histogram<TH1D> trapping_time_histo_;

        Histogram<TH1D> gain_histo_;

    };

} // namespace allpix