Module bca 

Amorphous binary-collision-approximation (BCA) transport in a layered or voxelised target, with full recoil cascades.

Model

The target is any Geometry: a Stack of layers or a VoxelGrid, whose regions are homogeneous, structureless (amorphous). A particle alternates straight free flights and binary elastic collisions with target atoms, losing energy continuously to electrons along each flight. This is the amorphous-target BCA of J. P. Biersack and L. G. Haggmark, Nucl. Instrum. Methods 174 (1980) 257 (TRIM), within the general BCA framework of M. T. Robinson and I. M. Torrens, Phys. Rev. B 9 (1974) 5008 and the treatment in W. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer, 1991). Everything here is implemented from those publications; no code from the Tier B/C programs named in CONTRIBUTING.md was consulted.

One history:

1. The primary enters at the geometry’s entry point (the origin of the front face of a stack; the centre of the front face of a voxel grid unless Bca::with_entry_point chooses another) with the beam’s direction.
2. Free flight of length tau * lambda (MeanFreePath), with the electronic loss N S_e(E) s taken at the energy at the start of the segment (Bragg additivity over the layer composition, crate::ion::stopping::bragg).
3. Collision with a partner drawn by stoichiometry, at impact parameter p = p_max sqrt(R) (uniform over the disc of radius p_max) and a uniform azimuth. The centre-of-mass angle comes from the ScatteringTable; lab angles and the transfer T from kinematics.
4. Recoil: if T > E_d (the partner’s displacement energy) the atom is displaced with energy T - E_b (the lattice binding E_b stays in the lattice). Otherwise T stays in the lattice at the site. This is the displacement criterion of Biersack and Haggmark (1980) and Eckstein (1991). Displaced atoms are followed in turn as full cascades (BcaConfig::follow_recoils), from an explicit stack (no recursion), and draw from the same per-history random stream, so the stream is keyed on the primary index only and results do not depend on the thread count.
5. A particle stops when its energy falls below its cutoff (BcaConfig::primary_cutoff_ev, BcaConfig::recoil_cutoff_ev), or leaves through the front or back face if it can overcome the planar surface barrier (kinematics::refract_out); otherwise it is reflected specularly back into the target.

Free-path convention and layer boundaries

Particle::layer is the particle’s region in the geometry: the layer index in a stack, the flat voxel index in a voxel grid. Material data (scattering, stopping, energies) is cached per material, shared by all voxels of that material. The surface barrier of an escaping atom uses the E_s of the material it leaves from and the actual outward normal of the face (see kinematics::refract_out_normal); a particle that cannot escape is reflected specularly about that normal. Escapes through the lateral faces of a voxel grid are reported as Face::Side and counted in EnergyBudget::lateral. crate::tally::ion::IonTally and the CLI tallies describe stacks only; use crate::tally::voxel::VoxelTally, SummaryTally or a custom BcaTally for voxel runs.

A flight is drawn as a dimensionless number of mean free paths tau and converted to a length with the local lambda. When a flight reaches a change of material (or a periodic wrap of a voxel grid) it is truncated exactly there, the electronic loss for the truncated length is applied, the material is switched, and the flight continues in the new region with the unused part tau - s/lambda (the flight is not redrawn). A voxel face between voxels of one material is not an event; the region a flight ends in is still tracked (Geometry::flight), so the collision, its tallies and its recoils carry the region of the collision site. The collision happens in the material where the flight ends, with a partner and p_max from that material. For exponential paths this is exact by memorylessness; for the constant path it means a layer split of one material changes nothing but floating-point rounding, which the split_layer_is_equivalent test checks.

Weak collisions (optional)

With the constant free path every flight ends in one collision with p <= p_max, and the nuclear loss of collisions with p > p_max is dropped, while the electronic loss of the flight is charged in full. At cascade-tail energies (a few to tens of eV) the dropped part is large: for Cu on Cu with the ZBL potential the disc p <= p_max carries 45 % of the nuclear stopping cross section at 5 eV and 54 % at 10 eV, so the electronic share of a cascade comes out far too high (issue #64; the level-1 check damage.cascade_electronic_share measures it).

BcaConfig::weak_collisions = K adds the “weak” collisions of W. Moller and W. Eckstein, TRIDYN - Binary collision simulation of atomic collisions and dynamic composition changes in solids, report IPP 9/64, Max-Planck-Institut fur Plasmaphysik, Garching (1988) (the long write-up of Comput. Phys. Commun. 51 (1988) 355), read at https://pure.mpg.de/rest/items/item_2131703/component/file_2131702/content. Read only the report body, PDF pp. 1-47: never Appendix 1 (the TRIDYN program listing, PDF pp. 48-86), which is Tier C (see CONTRIBUTING.md). Quoted:

• p. 14: “‘weak’ collisions might occur with more distant atoms which might contribute to energy loss and angular deflection. The present version allows up to three additional weak collisions with impact parameters larger than p_max, each of them representing one additional atomic volume”, with eq. (9) replaced by eq. (26), p_k^weak = p_max sqrt(k + r_p), k = 1, 2, 3 (r_p uniform in [0, 1)): the k-th partner is uniform over the annulus between p_max sqrt(k) and p_max sqrt(k + 1), of area pi p_max^2. Fig. 3 (p. 11) draws the first-order partner at its own azimuth.
• p. 26: “the weak collision loop is entered which finally defines new directions after each of the simultaneous collisions. (Actually, the last passage of the weak collision loop represents the hard collision.) A primary recoil may be generated and stored for each hard collision.”
• p. 35: the program sums the elastic energy transfers of the weak and hard collisions, and their local inelastic (Oen-Robinson) losses.
• Figs. 7 and 8, pp. 29-30 (flow charts of the projectile and cascade loops): inside the weak-collision loop, after the partner’s species, azimuth and impact parameter are drawn, a “Target beyond Surf.?” test skips the passage when the partner would lie outside the target; and “Elastic Energy Loss” is applied once, after the loop.

So each collision step here is: the K weak collisions, k = 1..K, each with a partner drawn by stoichiometry, p = p_max sqrt(k + R) and a uniform azimuth, then the hard collision with p = p_max sqrt(R). Each deflects the particle and takes its transfer T; under ElectronicLoss::EquipartitionLsOr each also takes its Oen-Robinson local loss, so the local half is sampled out to p_max sqrt(K + 1). A weak collision never makes a recoil: its T stays in the lattice at the site (LatticeDeposit::Weak, counted in EnergyBudget::lattice). Random draws per weak collision: partner, R, azimuth, in that order and before the hard collision’s, so K = 0 draws exactly what the engine drew before the option existed and reproduces it bit for bit.

Surface test. As in Figs. 7 and 8, a weak collision whose partner would lie in front of the front surface (in vacuum) is skipped; its draws are still made. The partner sits at distance p from the path on the side opposite to the deflection (where a hard collision sends its recoil). In TRIDYN the test also covers the hard collision (the loop’s last passage); here it does not, because that would change the engine without weak collisions, which this option leaves untouched. Whether the hard collision should be tested too is left as a separate model decision. The back face of a finite target is not tested: the flow charts test “the surface” only.

Deviation from TRIDYN: the energy of each collision. The flow charts apply the elastic loss once, after the loop, and p. 35 sums the transfers, so in TRIDYN every collision of a step is evaluated at the energy the step starts with. Here each is evaluated at the energy left after the previous one (and in the direction after it), so no transfer can exceed the particle’s energy. Evaluated at the starting energy, the summed transfers of a slow atom often exceed its energy, and the report shows no rule for that case: for Ar 1 keV on Cu with E_d = E_s and K = 3, 18 % of the binary collisions would need a clamp, about 200 eV per ion in all, and the sputter yield would depend on how the clamp is applied. Measured on that problem (20 000 ions, seed 1): 0.882 as here; 0.875 with only the weak collisions at the starting energy (each transfer clamped to the energy left, the hard collision at the energy left); 0.976 with every collision at the starting energy and each transfer clamped in loop order. Since the clamp is not in the source, the sequential choice is kept.

With K weak collisions the nuclear loss per path is the stopping cross section integrated to p_max sqrt(K + 1) instead of p_max (93 % of the full value at 5 eV for Cu on Cu with K = 3). It is available with MeanFreePath::Constant only (p_max is the constant-path radius); MeanFreePath::EnergyDependent rejects it, since that convention drops small-angle collisions by design.

Default 0, because there is no single published convention: Biersack and Haggmark (1980), whose constant free path this engine follows, use one collision per flight; TRIDYN allows up to three (the report read does not give a default); SDTrimSP, a TRIDYN descendant, defaults to two for projectiles and two for recoils (iwc, iwcr: A. Mutzke et al., SDTrimSP Version 6.00, report IPP 2019-02 (2019), Table 13, p. 73, “number of ring cylinders for weak simultaneous collisions”); and the RustBCA manual documents weak_collision_order with a default of 0 (its GitHub wiki, “Standalone Code: Input File”; manual only, Tier B).

The option is not a free improvement, which is the other reason it is off by default. It brings the electronic share of a cascade in line with the Lindhard partition integral equation solved for the engine’s own inputs (level-1 rows damage.cascade_electronic_share.*: for Cu on Cu at 1 keV about 0.79 without weak collisions and 0.35 with K = 3, against 0.38). But a weak transfer stays in the lattice, and at a few eV it is a large part of the energy (Cu on Cu at 5 eV, ZBL: on average 46, 25 and 13 % in the first three annuli), so slow atoms near the surface stop sooner: the Ar 1 keV on Cu sputter yield with E_d = E_s falls from 1.93 to 0.88 with K = 3 (0.66 without the surface test; docs/validation.md, level 2). Use it for the energy partition (damage energy, electronic share); its effect on yields is large and has not been checked against measured yields.

Energy bookkeeping

Every energy change is subtracted from the particle and added to exactly one field of an EnergyBudget, so each history conserves energy up to rounding (EnergyBudget::residual). Electronic losses larger than the particle’s energy are clamped to it.

Vacancy, interstitial and replacement counts, and the NRT damage estimate, are kept by the tally crate::tally::IonTally from the events reported here (crate::ion::damage has the models).

Crystal regions

Everything above is the amorphous partner model. Bca::with_crystal switches chosen regions to the crystal flight model (partners from an explicit lattice, in crystal); the other regions, and engines that do not call it, are untouched and bit-identical to the amorphous engine.

Not modelled here (see the issue tracker)

Damage accumulation in crystal regions (the lattice there is perfect, static or thermally vibrating; the local Oen-Robinson loss of ElectronicLoss::EquipartitionLsOr is taken at every lattice partner, see crystal), target composition changes with fluence, and refraction of the incident beam at the entrance surface (negligible at keV energies).

Re-exports

pub use crystal::CrystalMetadata;

pub use crystal::CrystalTarget;

pub use crystal::Thermal;

pub use crystal::ThermalMetadata;

pub use tally::BcaTally;

pub use tally::ElectronicChannel;

pub use tally::EnergyBudget;

pub use tally::Face;

pub use tally::LatticeDeposit;

pub use tally::SummaryTally;

Modules

crystal

The crystal flight model: collision partners taken from explicit lattice sites instead of a random impact parameter (step 20b of the M2 crystal plan), optionally displaced by thermal vibration (step 21b).

kinematics

Two-body elastic kinematics, direction rotation and planar surface refraction. Pure functions, no state.

tally

What the BCA engine reports, and a small built-in summary tally.

Structs

Bca

The BCA engine for one beam, target and configuration.

BcaConfig

Run configuration.

Beam

The incident beam.

DivergenceMetadata

Run metadata of a beam divergence (Bca::divergence_metadata): the distribution, its width and the incidence policy, so a result states the spread it was run with.

HistoryBuffers

Reusable working memory for Bca::history_in: the stack of particles still to be followed and the per-flight scratch. A history that runs in buffers which have already seen a cascade as large allocates nothing, so a worker thread allocates only while its buffers grow to the largest cascade it meets, not per history and not per collision.

Particle

A moving or stopped particle.

Enums

BcaError

Errors from setting up or running the engine.

ElectronicLoss

Electronic energy-loss treatment.

HistoryError

Why a single history failed.

MeanFreePath

How free-flight lengths and the impact-parameter limit are chosen.

Constants

DIVERGENCE_STREAM_WORD

Word position (in the 32-bit words of the ChaCha stream) of the segment of a history’s random stream that beam-divergence draws come from: 2^65. The transport draws start at word 0 and never reach it, and it is disjoint from the thermal-displacement segment (2^66) and the crystal-shift segment (2^67) (each is far shorter than the gap to the next), so enabling divergence leaves every other draw of the history unchanged.

MAX_DIVERGENCE_ATTEMPTS

Most deflected directions drawn for one primary before Bca::with_divergence’s inward conditioning gives up with BcaError::BeamDivergence. For the small-angle spreads the input front end accepts, one draw is inward with probability about one half or more, so exhausting 1000 attempts has probability below 2^-1000; the bound only guarantees that a pathological setting cannot hang a run.

MAX_WEAK_COLLISIONS

Largest BcaConfig::weak_collisions: “up to three additional weak collisions” (Moller and Eckstein, IPP 9/64 (1988), p. 14).