# RREA-Feedback-Geant4 Geant4 simulation and analysis for estimating the relativistic runaway electron avalanche (RREA) multiplication *M* and the relativistic-feedback factor *γ* in a realistic, altitude-dependent thundercloud electric field and atmosphere, and for locating the self-sustaining feedback threshold *γ* = 1. Direct transport of a full avalanche becomes infeasible above threshold (*M* grows like *e^ξ*, with *ξ* the number of avalanche lengths). The analysis therefore uses a factorized estimator: it measures the per-electron ingredients (avalanche length, carrier yield, seed effectiveness, response to single feedback carriers) in short runs and assembles *M*(*E₀*) and *γ*(*E₀*) on a grid of peak fields *E₀*. ## Contents ``` CMakeLists.txt, src/ Geant4 application rrea_feedback (C++, NRLMSISE-00 in C) scripts/ build/run wrappers, calibration runs, factorized analysis, self-tests data/ field profile and charge-layer model docs/factorized/ calibration and result files (see "Results") CITATION.cff, LICENSE ``` ## Simulation setup * **Electric field.** A vertical profile −*E_z*(*z*) (`data/efield_profile.csv`, generated by `scripts/generate_efield_profile.py` from the charge-layer model `data/lma_charges_tweaked.json`; unscaled peak 125.54 kV/m), rescaled so that its maximum equals `--e0-kv-m`. Its column `minus_Ez_kV_per_m` is −*E_z*: positive values are a downward field, which accelerates electrons upward. The field fills a vertical cylinder (`--field-radius-km`, default 5 km). `data/lma_charges_adjusted_peak121p3.json` and `data/efield_profile_peak121p3.csv` are the same model with the upper positive layer moved so that the profile peaks at 121.3 kV/m. * **Atmosphere.** Layers of 50 m with NRLMSISE-00 densities at 28.7° N, −80.85° E, 29 July 2023, 20:30:30.640 UTC, F10.7 = F10.7A = 160, Ap = 7. * **Seeds.** 1 MeV electrons moving upward, sampled uniformly over the cylinder cross-section just above the lower boundary of the region where the field exceeds the RREA threshold (`--seed-altitude-mode lower-boundary`). * **Physics.** `G4EmStandardPhysics_option1` with its default stepping, and `G4DecayPhysics`. Production cuts are set with `--production-cut-mm` and `--efield-production-cut-mm`. * **Feedback.** Returned photons and positrons that would seed new avalanches are recorded; feedback is truncated after one generation (`--max-feedback-generation 1`), so a run measures one-generation estimators. `build/rrea_feedback --help` lists all options. ## Requirements * Linux. * Geant4 10.7 (tested with 10.7.4), CMake ≥ 3.16, a C and a C++17 compiler. * Python 3 (standard library; matplotlib is optional and only used for plots). ## Build and run ```bash scripts/build_local.sh # uses geant4-config from PATH, or G4_PREFIX= scripts/run_local.sh --seeds 100 --e0-kv-m 120 --field-radius-km 5 --output-dir output/e0_120 ``` `run_local.sh` sources `geant4.sh` from the same installation, clears data-set variables left over from other Geant4 installations, and builds the executable if needed. Each run writes, in its output directory: * `feedback_summary.txt` / `feedback_summary.csv`: run settings and results, including `gamma_effective` (the one-generation feedback factor: upper-plane electrons produced by generation-1 feedback seeds divided by those produced by the original seeds) with event-bootstrap 95% intervals; * `plane_crossings.csv`, `event_plane_crossings.csv`, `event_stats.csv`: upward electron crossings of altitude planes, per run and per event; * `rrea_threshold_profile.csv`: density, field, RREA threshold and friction versus altitude. ## Factorized analysis The full calibration and analysis chain runs with ```bash bash scripts/run_recalibration_campaign.sh # MAXJOBS=8 SEEDS_SLAB=120 ... to override ``` in four phases: P, one-seed profile runs on the *E₀* grid; C1, slab calibration of the avalanche length and carrier yield (`run_slab_calibration.sh`, `fit_lambda_model.py`, `build_carrier_yield.py`) with the same-field cross-altitude check `check_altitude_gate.py`; C2, 1 MeV seed-effectiveness anchors (`factorized_gamma.py build-eps`) with the check `validate_gates.py`; C3, carrier sampling, response runs and assembly (`factorized_gamma.py sample`, `run_response.sh`, `factorized_gamma.py assemble`). It skips runs that are already complete. Each Python script lists its options with `--help`. `scripts/build_eps_direct.py` calibrates the seed effectiveness directly against transported multiplication. ```bash python3 scripts/selftest_factorized.py # fast self-tests of the analysis code ``` ## Results * `docs/factorized/lambda_g4.json`: Geant4 avalanche length versus reduced field and density, from slabs at 11 field values between 86 and 270 kV/m (fit `λ = C1/(E − C2 ρ̂)` with C1 = 9411 kV and C2 = 246 kV/m). The slab field profiles are regenerated in `data/slabs/` by `run_slab_calibration.sh`. * `docs/factorized/carrier_yield.json`: density-reduced feedback-carrier yield. * `docs/factorized/eps_seed.json`: 1 MeV seed effectiveness (asymptote 2.44, 95% interval 1.37 to 3.65). * `docs/factorized/factorized_results_mu1_300k.csv`: *M*_est and *γ*_eff with 95% intervals at *E₀* = 110, 120 and 130 kV/m (all photon directions, 3 × 10⁵ carriers per point). The retained feedback grid combines the 110 and 130 kV/m points of that file with runs of 10⁶ carriers at 119 and 120 kV/m (all photon directions): | *E₀* (kV/m) | carriers | *γ*_eff [95% interval] | photon | positron | |---:|---:|---:|---:|---:| | 110 | 3 × 10⁵ | 0.125 [0.066, 0.212] | 0.063 | 0.062 | | 119 | 10⁶ | 1.301 [0.699, 2.015] | 0.808 | 0.493 | | 120 | 10⁶ | 1.399 [0.769, 2.201] | 0.732 | 0.667 | | 130 | 3 × 10⁵ | 11.26 [5.84, 18.21] | 5.77 | 5.50 | Log-linear interpolation gives *γ* = 1 at *E₀* ≈ 118.0 kV/m, and 116.2 to 121.3 kV/m from the upper and lower 95% bounds. These files were produced by this code before its physics list was set to `G4EmStandardPhysics_option1` (it then registered `G4EmStandardPhysics_option4`); they were not regenerated. The raw Geant4 run outputs behind them are not included. ## Verification With Geant4 10.7.4 the code builds without warnings, the self-tests pass, and a one-seed run at *E₀* = 85 kV/m completes. ## Not included Raw run outputs (including the response runs behind the table above), other seed energies, the earlier direct-transport study and its tools, and the method documents and report of the original project. Some code comments still refer to those documents. ## Third-party code NRLMSISE-00 C version by Dominik Brodowski (`src/src/msis/C/`), from the model of M. Picone, A. Hedin and D. Drob. ## License and citation MIT (see `LICENSE`). See `CITATION.cff`. ## Contact David Sarria, Department of Physics and Technology, University of Bergen, Norway (david.sarria@uib.no).