pAvalancheMC: RPC Monte Carlo Simulation

• pAvalancheMC is a parallel Monte Carlo avalanche simulation tool integrated in Garfield++ for RPC studies with dynamic space‐charge modifications.
• It uses OpenMP parallelization and a line-charge representation to efficiently model electron transport, multiplication, and induced signals.
• The method combines image-charge corrections, anisotropic diffusion, and realistic field updates to enhance detector performance simulation.

Searching arXiv for the specified pAvalancheMC papers and closely related RPC simulation work. I’ll gather the relevant arXiv records for the Garfield++/RPC pAvalancheMC work and its thesis follow-up before writing the article. Searching arXiv now. pAvalancheMC is a parallel Monte Carlo avalanche class introduced in Garfield++ for Resistive Plate Chamber (RPC) simulation in regimes where dynamic space-charge effects materially modify the applied field. It was developed to make 3D particle-based avalanche calculations tractable in operating domains that include avalanche growth, saturated avalanche, and streamer-like onset, while remaining integrated with standard Garfield++ components such as neBEM, MAGBOLTZ, TrackHeed or HEED, and ROOT random-number infrastructure. In its published form, pAvalancheMC combines OpenMP parallelization with a line-charge representation of the evolving charge cloud; in the later thesis treatment, this formulation is extended by image-charge modeling for layered dielectric electrodes, so that the total field can include applied, space-charge, and polarization contributions (Dey et al., 2022, Dey, 29 Jul 2025).

1. Origin, scope, and detector context

pAvalancheMC was introduced to address a specific deficiency in conventional Garfield++ avalanche tools for RPC studies: large avalanches are computationally expensive, while the evolving electron and ion cloud can substantially distort the local electric field and thereby alter transport, multiplication, and signal formation. The class targets RPC operating regimes in which the number of charges becomes large enough to modify the applied field, including the transition from ordinary avalanche growth to saturated avalanche and to streamer-like behavior (Dey et al., 2022).

The detector context is the parallel-plate RPC, including single-gap and timing RPC variants with resistive or dielectric electrodes such as bakelite or glass, graphite coating, and segmented pickup strips. The thesis formulation explicitly situates the method in the ICAL/INO program, where glass RPCs operate predominantly in avalanche mode and detector optimization depends on accurate estimates of timing, efficiency, charge spectra, and rate capability under space-charge and electrode-polarization effects (Dey, 29 Jul 2025).

A central distinction of pAvalancheMC is that it remains a Monte Carlo transport-and-multiplication engine rather than a fully microscopic collision-by-collision simulation. This allows substantially faster event generation than AvalancheMicroscopic while retaining field-dependent transport, stochastic multiplication, and induced-signal calculation. The intended application domain is therefore detector-level performance simulation in parameter regimes where self-consistent field updates are indispensable.

2. Field model, transport physics, and avalanche growth

The basic field decomposition in the original Garfield++ implementation is

$E_t = E_a + E_s,$

where $E_a$ is the applied field and $E_s$ is the dynamically recomputed space-charge field. In the thesis extension, the total field includes an additional polarization term from image charges in the dielectric stack, so that transport can be updated under the modified local field at every time step (Dey et al., 2022, Dey, 29 Jul 2025).

Transport and multiplication coefficients are taken from MAGBOLTZ. The quantities used are the drift velocity $v_d(E)$, longitudinal and transverse diffusion coefficients $D_L(E)$ and $D_T(E)$, Townsend ionization coefficient $\alpha(E)$, and attachment coefficient $\eta(E)$. The diffusion model is Gaussian and anisotropic under field. In one formulation,

$$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$

while the thesis also gives the anisotropic Gaussian forms

$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$

and

$E_a$0

Avalanche growth is governed by the effective Townsend coefficient $E_a$1. The expected mean electron population along a drift path is written as

$E_a$2

and the stochastic implementation uses a modified Yule–Furry process with attachment. In the thesis description, each drift step is subdivided so that the subdivision length is approximately $E_a$3, with probabilities $E_a$4 and $E_a$5 for ionization and attachment, respectively (Dey et al., 2022, Dey, 29 Jul 2025).

This formulation makes saturation a direct field-feedback effect rather than an imposed cap. As space charge accumulates, $E_a$6 is reduced in parts of the avalanche, lowering $E_a$7 and increasing the relative importance of $E_a$8. When $E_a$9 becomes negative in most voxels, net growth stalls. In Ar/CO$E_s$0 (80/20), the reported saturated avalanche shows a knee around $E_s$1 ns at an electron count of approximately $E_s$2, close to the Raether limit traditionally associated with streamer transition (Dey et al., 2022).

3. Software architecture and algorithmic workflow

pAvalancheMC is built as a Garfield++ class that interoperates with the standard RPC simulation stack. The component roles described in the sources are summarized below.

Component	Function	Notes
neBEM	Static and weighting fields	OpenMP-parallelized; FastVol interpolation
MAGBOLTZ	Transport and multiplication tables	$E_s$3, $E_s$4, $E_s$5, $E_s$6, $E_s$7
TrackHeed / HEED	Primary ionization clusters	Positions and generation times
ROOT TRandom3	Per-thread random streams	Mersenne Twister, period $E_s$8
OpenMP	Parallel execution	Drift, gain, field update, signal accumulation

A typical run initializes geometry and static fields, loads transport tables, configures OpenMP threads and per-thread TRandom3 instances through SetNumberOfThreads(n), generates primary clusters, and acquires weighting-field maps. Primaries are sorted by generation time and started when the simulation clock satisfies $E_s$9. The time-stepping loop then advances electrons, applies diffusion, samples multiplication and attachment, updates the charge lists, and recomputes the field when the charge threshold for space-charge activation is exceeded (Dey et al., 2022).

The configuration described in the published examples includes a time step $v_d(E)$0 ps, a space-charge activation threshold around $v_d(E)$1, and cylindrical segmentation parameters with typical values $v_d(E)$2 cm and $v_d(E)$3. Garfield++ interfaces named in the sources include SetTimeStep(20e−12), EnableSpaceCharge(true) or SetSpaceChargeOn(true), SetSpaceChargeThreshold(1e4), SetVoxelization(dz, dr, dphi) or SetGrid(zStep, rStep, phiStepDeg), SetImageChargeOn(true/false), SetImageCut(percent), SetElectrodeProperties(...), and RunAvalanche(...) (Dey et al., 2022, Dey, 29 Jul 2025).

OpenMP parallel regions are applied to electron drift and diffusion, multiplication and attachment sampling, space-charge field summation, and induced charge or current accumulation. The per-thread random-number strategy is explicit: TRandom3 instances are allocated per thread with distinct seeds to avoid correlated random streams.

4. Space-charge representation and dielectric polarization

The distinctive numerical device in pAvalancheMC is the cylindrical voxelization of the avalanche charge cloud. Electrons and ions are binned into $v_d(E)$4 slices in $v_d(E)$5, $v_d(E)$6 concentric radial rings, and $v_d(E)$7 azimuthal sectors with $v_d(E)$8. Each occupied voxel is represented by a straight line-charge segment with constant linear charge density $v_d(E)$9, and curved azimuthal segments are approximated by straight segments of length $D_L(E)$0 (Dey et al., 2022).

In the original implementation, the field of a finite line segment is evaluated numerically through

$D_L(E)$1

and voxelwise field components $D_L(E)$2, $D_L(E)$3, and $D_L(E)$4 are accumulated. The paper notes that $D_L(E)$5 is generally small, but not identically zero; one example gives a maximum $D_L(E)$6 of $D_L(E)$7 at $D_L(E)$8 ns, indicating axial asymmetry (Dey et al., 2022).

The thesis extends this by introducing a segmented-line model together with image charges in layered dielectrics. The RPC stack is treated as a five-layer system: graphite, dielectric electrode, gas, dielectric electrode, graphite. Reflection and equivalence factors at interfaces are

$D_L(E)$9

with $D_T(E)$0 for a conductor. Image charges are generated iteratively and truncated using a percentage-change criterion,

$D_T(E)$1

with a representative cut of $D_T(E)$2 (Dey, 29 Jul 2025).

The same source gives the equivalent RPC relaxation time

$D_T(E)$3

and reports the example $D_T(E)$4 mm, $D_T(E)$5, $D_T(E)$6, yielding $D_T(E)$7 s. The static polarization field is thus incorporated through the image series, while the full time-domain polarization evolution $D_T(E)$8 is not implemented (Dey, 29 Jul 2025).

A practical implication is that pAvalancheMC evolved from a solver for $D_T(E)$9 to a more general RPC field engine in which $\alpha(E)$0 can include dielectric-image contributions. This suggests that later usage is not limited to gas-phase self-fields alone, but also addresses the layered electrostatics of realistic RPC stacks.

5. Signal induction, validation, and numerical performance

Induced signals are computed through the Ramo–Shockley formalism using weighting fields from neBEM. The induced charge is written as

$\alpha(E)$1

and the induced current, when requested, as

$\alpha(E)$2

The thesis also gives an analytic approximation for the weighting field in a single-gap three-layer RPC geometry,

$\alpha(E)$3

These signal calculations are integral to timing RPC studies rather than an external post-processing stage (Dey et al., 2022, Dey, 29 Jul 2025).

A timing RPC example with a $\alpha(E)$4 mm gap at $\alpha(E)$5 kV/cm reports induced-charge distributions at $\alpha(E)$6 V, $\alpha(E)$7 V, and $\alpha(E)$8 V fitted with a Polya function,

$\alpha(E)$9

The mean charge parameter increases from approximately $\eta(E)$0 fC to $\eta(E)$1 fC to $\eta(E)$2 fC as the voltage increases, and the parameter $\eta(E)$3 shifts from $\eta(E)$4 to $\eta(E)$5 to $\eta(E)$6. The low-charge inefficiency peak was excluded from the fits, and the trend is reported as consistent with experimental observations attributed to Fonte and Lippmann (Dey et al., 2022).

The space-charge dynamics were validated against earlier solvers and against experimental trends. The 3D Cartesian implementation reproduces features reported in earlier 2D rotationally symmetric models, including radial blow-up and contraction and comparable gain dynamics. The thesis additionally states that the segmented-line field matches the ring model after a $\eta(E)$7 scaling and that source-plus-image fields agree with neBEM near the electrode when image charges are included, with ratios approximately unity there and about $\eta(E)$8 to $\eta(E)$9 elsewhere where the field is not vanishing (Dey, 29 Jul 2025).

Parallel performance is a defining feature. In a timing RPC of area $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$0 cm$$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$1, gas gap $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$2 mm, electrode thickness $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$3 mm, and field $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$4 kV/cm, the speed-up $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$5 for $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$6 avalanches was fitted as

$$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$7

with asymptotic speed-ups $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$8 without space charge and $$\sigma_L = D_L \sqrt{v_d \delta t}, \qquad \sigma_T = D_T \sqrt{v_d \delta t},$$9 with space charge. The saturation of speed-up is attributed to OpenMP overheads and the serial fraction of the code. neBEM parallelization was tested up to $$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$0 cores while preserving precision and significantly reducing field computation times (Dey et al., 2022).

6. Assumptions, limitations, and research outlook

The modeling assumptions are explicit. Space charge is represented by straight line segments on a cylindrical grid; the field is treated as uniform within each voxel; near-field singularities are moderated by finite-segment integrals and voxel averaging; and positive ions and negative ions are static unless ion mobility is supplied. This last simplification is justified in the sources by the large disparity between electron drift velocity and ion mobility on the simulated time scales (Dey et al., 2022).

The code does not implement full streamer physics. It can identify streamer-like onset through rapid gain divergence and high local fields, but photon transport and photoionization are absent. This is important because saturated avalanche and streamer onset are treated differently: the former corresponds to field suppression such that $$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$1 over most of the avalanche volume, whereas the latter is inferred when electron number and $$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$2 rise sharply without a clear saturation plateau. In C$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$3H$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$4F$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$5/i-C$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$6H$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$7/SF$$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$8 at $$\phi_L(z,l) = \frac{1}{(\sqrt{2\pi l}) D_L} \exp\!\left(-\frac{(z-z_0)^2}{2 D_L^2 l}\right),$$9 kV/cm, the reported regime is characterized by a peak followed by a sudden increase in electron number while $E_a$00 exceeds approximately $E_a$01–$E_a$02 kV/cm; full streamer simulation is deferred to future work (Dey et al., 2022).

The thesis adds further caveats: ion–ion interactions, detailed electrode surface roughness, graphite resistivity nonuniformity, and a complete transient polarization model are simplified or absent. Sensitivity to electrode relative permittivity $E_a$03, thickness $E_a$04, and gas gap $E_a$05 is emphasized: the field increases rapidly with $E_a$06 and saturates beyond about $E_a$07 for fixed $E_a$08 and $E_a$09, decreases with increasing $E_a$10 for fixed $E_a$11 and $E_a$12, and decreases with increasing $E_a$13 for fixed $E_a$14 and $E_a$15 (Dey, 29 Jul 2025).

Within these bounds, pAvalancheMC functions as a detector-physics tool for RPC optimization rather than a universal plasma-discharge simulator. Its established strengths are self-consistent avalanche growth under dynamically updated fields, induced-signal prediction using realistic weighting fields, and computational tractability through OpenMP parallelism. The stated research direction is toward photon transport, more detailed time-dependent electrode charging, richer diffusion and field solvers, and broader acceleration strategies such as GPU or MPI-based parallelization (Dey, 29 Jul 2025).