Module 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).

Choosing the model per region

Geometry is unchanged. A Bca is amorphous everywhere unless Bca::with_crystal names a CrystalTarget (lattice and orientation) and the list of regions (layers of a stack, voxels of a grid) it fills. A particle whose current region is in that list uses the crystal model; in every other region, and for a Bca without with_crystal, the amorphous model is used unchanged (and bit for bit: the amorphous code path draws the same random numbers and does the same arithmetic as before; the test tests/crystal_amorphous_identity.rs compares complete reports). Several calls give several crystals (different lattices or orientations) in different regions.

The material of a crystal region still supplies the displacement and binding energies, the stopping and the surface barrier; it must contain every element of the lattice and its atom density must agree with the lattice’s to 5 %. The crystal model needs the constant-path convention and no weak collisions (the lattice replaces both). Beam direction and surface stay in the lab frame of the geometry; the orientation’s lab-to-crystal rotation (Orientation::to_crystal) is the only thing read from it (the beam is the Beam’s, so give the Beam the orientation’s tilt and twist, as Orientation::beam_lab does).

Lattice frame of a history

The crystal position of a lab point r is R r + u, with R the orientation’s rotation and u a translation drawn uniformly from the rectangular cell that tiles the crystal, [0, e_x) x [0, e_y) x [0, e_z) with the edges of Lattice::orthogonal_cell (the cube [0, a)^3 for a cubic lattice), once per primary history (three uniform numbers per crystal, drawn from a copy of the history’s stream positioned at word 2^67, so the transport draws are the same with or without crystals), shared by the primary and all its recoils. A beam of finite width samples all positions of the unit cell, so this is the static equivalent of a random entry point; it also makes the surface termination random. Lattice sites that lie outside the target (in front of the surface) or in a region that is not part of this crystal are not partners.

Collision sequence (after DISPLATH, Tier A)

The ion moves in straight segments. For a segment from position r0 along the unit direction d the partners are the lattice sites within p_max of the line, ordered by the path distance s = (r - r0) . d of their point of closest approach (crate::ion::crystal::LatticeSearch). The scheme is the one of DISPLATH (https://github.com/permissionx/DISPLATH, MIT, commit 7f461141c518305e6a5de8ce7e37dd314cd20ccb, src/dynamics.jl: GetTargetsFromNeighbor and Collision!, and SimultaneousCriteria in src/geometry.jl), re-implemented here on LatticeSearch with the lattice-site list of this crate, and with the project’s own scattering tables, kinematics and energy bookkeeping. DISPLATH is Copyright (c) 2024 裴幂許Permission (MIT); the notice is in THIRD_PARTY_LICENSES.md.

• Free path. The ion flies straight to the closest-approach point of the nearest partner (smallest s, ties broken by the search’s documented order). There is no free-path length parameter: the path is the crystal’s. The nonlocal electronic loss of the model chosen in the BcaConfig is taken over that length at the energy at the start of the segment, as in the amorphous model. Under ElectronicLoss::EquipartitionLsOr that is half of the Lindhard-Scharff loss, and the other half is local (next bullets and “Local electronic loss”). When no partner lies within search_length the ion flies that distance and searches again (a channel). The geometry is asked about every segment, so surfaces and interfaces cut it exactly like an amorphous flight.
• Simultaneous collisions. Further sites join the nearest one if, for every partner t already accepted, Δ = s - s_t (≥ 0) satisfies Δ <= q_max, p_t² + Δ² <= p_max² and p² + Δ² <= p_max² (DISPLATH’s SimultaneousCriteria: sites at similar path distance and inside a sphere of radius p_max around the collision point). DISPLATH sets q_max to the sum of the two atomic radii (about the nearest-neighbour distance); here the default is a small fraction of it, see the next section.
• Collision. Each partner i is scattered from the incoming state at its own impact parameter p_i (the distance from the site to the line) with the scattering table: transfer T_i = γ E sin²(Θ_i/2), the recoil leaving along sin(Θ_i/2) d + cos(Θ_i/2) n_i, with n_i the unit vector from the line to the site (the half-angle form of kinematics::lab_angles: the recoil angle is (π - Θ)/2). The projectile takes the momentum that is left, P = P0 - Σ sqrt(2 m_i T_i) r_i, and the energies are scaled by the common factor λ = (E - Σ Q_i) / (|P|²/2m + Σ T_i) so that energy is conserved with the local electronic losses Q_i (DISPLATH’s Collision_!). For one partner and no local loss this is the ordinary binary collision (λ = 1). The projectile stays at the closest-approach point of the nearest partner; the time-delay offsets of Robinson and Torrens (1974) are not applied (DISPLATH does not apply them either).
• Recoils start at their lattice site, with the usual displacement criterion (BcaConfig::follow_recoils) and follow the same model. The sites that were just hit are excluded from the next search (so is the site the recoil starts on): the ion cannot collide twice in a row with the same atom.

Choice of the search parameters

Measured on this engine (B 5 keV into Si, 4000 ions, ZBL, Lindhard-Scharff, 7 degrees tilt, 22 degrees twist, static lattice; ratios of Rp to the amorphous result), with nn the nearest-neighbour distance (2.35 Å):

• With q_max = nn (DISPLATH’s rule) the result depends strongly on p_max: p_max = 1.53 Å gives Rp ratio 1.11, but p_max = 2.3 Å gives 4.4. The merged collisions of distant sites cancel the kicks of a channel’s walls, so almost every ion channels, at 20 and 30 degrees tilt as well. That is an artefact of merging over a sphere as large as the channel, not physics.
• With q_max <= 0.1 nn the result is independent of q_max and, for p_max >= 1.5 * 1.53 Å = nn, of p_max within the statistics (Rp ratio 1.08-1.09 for p_max = 2.3, 3.1 Å and q_max = 0.01-0.1 nn; at q_max = 0.2 nn the artefact starts).

So the defaults are p_max = nn and q_max = 0.05 nn: converged in both, and the cost grows with p_max. These are properties of this engine’s implementation of the criterion, checked on silicon only; other lattices need the same convergence check.

Measured checks (static lattice, 4000 ions, seed 1, Si, ZBL, Lindhard-Scharff)

tests/crystal_bca.rs (the statistical ones are #[ignore]d; release build). Ratios are crystal over amorphous:

Case	Rp	dRp	Notes
B 5 keV, 30° tilt, 17° twist	0.97	1.09	within 10 % / 15 %
As 30 keV, 30° tilt, 17° twist	0.93	1.13	within 10 % / 15 %
B 5 keV, 7° tilt, 22° twist	1.06	1.24	median 23.8 vs 23.8 nm, 99th percentile 76.2 vs 56.8 nm
As 30 keV, 7° tilt, 22° twist	1.01	1.88	median 22.2 vs 25.2 nm, 99th percentile 112.2 vs 50.8 nm
B 5 keV along <110>	6.1 (121 vs 20 nm)	3.2	91 % of ions deeper than twice the amorphous Rp (amorphous 5 %)

The random-direction check (Rp within 10 %, dRp within 15 % of the amorphous result) is made at 30°/17° only. That beam is at least 8.0° from every {100}, {110} and {111} plane. It passes within 1.3° of a {211} and 0.9° of a {311} plane, which are higher-index planes with a lower atomic areal density. The 7°/22° rows are not a random-direction check; see “The 7°/22° criterion” below.

Measured checks with thermal vibration (same set-up, THETA_D_SI, zero-point term on)

tests/crystal_thermal.rs (statistical ones #[ignore]d; release build). “Tail” is the fraction of primaries deeper than twice the amorphous Rp (amorphous: 0.034 at 7°/22°, 0.041 at 30°/17°), with its binomial error.

Case	Static	0 K	300 K	600 K
B 5 keV along <110>: Rp, nm	121.0	111.7	94.3	75.8
B 5 keV along <110>: tail	0.910	0.876 ± 0.005	0.833 ± 0.006	0.758 ± 0.007
B 5 keV, 30°/17°: Rp / dRp ratio	0.97 / 1.09	0.99 / 1.15	1.00 / 1.18	0.99 / 1.15
As 30 keV, 30°/17°: Rp / dRp ratio	0.93 / 1.13		0.95 / 1.17
B 5 keV, 7°/22°: Rp / dRp ratio	1.06 / 1.24	1.13 / 1.41	1.14 / 1.43	1.10 / 1.29
As 30 keV, 7°/22°: Rp / dRp ratio	1.01 / 1.88		1.05 / 1.80

The static and 300 K entries of the 30°/17° and 7°/22° rows were measured again for #225 (aarch64 macOS, commit 8513cfe). The 0 K and 600 K entries and the <110> rows are from PR #224. The crystal path is not bit identical across platforms (docs/architecture.md). The re-measured values differ from those of PR #224 by at most 0.02 in a dRp ratio. That is within the seed-to-seed spread below.

Along <110> the channeled tail falls with temperature, each step by more than five binomial standard errors (the #181 criterion asks for a monotonic decrease beyond the statistics). At the random direction the Rp ratio stays within 10 % at 300 K.

Off-axis directions get a longer tail with vibration, not a shorter one. At 7°/22° (B) and 30°/17°, the tail beyond twice the amorphous Rp grows from the static value as the amplitude rises from zero. It peaks near u1 ≈ 0.045-0.065 Å (0-300 K for Si) and falls again at higher amplitudes. PR #224 measured, for B at 7°/22° with 4000 ions: tail 0.080 static, 0.075 at u1 ≈ 0.0004 Å, 0.108 at 0 K, 0.115 at 300 K, 0.094 at 600 K, and 0.076 at 1500 K. The dependence is smooth in the amplitude. It is converged in the search parameters at 300 K: q_max = 0.01, 0.05 and 0.1 nn and p_max = 1.5 nn spread the dRp ratio over 1.39-1.45 and change Rp by under 1 %. So it is a property of the model as built, not of the thermal code path.

#225 measured how stable this is. The rise is significant at 16000 ions (B 5 keV), under both electronic-loss treatments. The amorphous reference uses the same loss as the crystal run (off_axis_tail_grows_with_vibration_under_both_losses):

B 5 keV, 16000 ions	Amorphous tail	Static tail	300 K tail	Rise	300 K Rp / dRp ratio
7°/22°, NonLocal	0.0350	0.0806 ± 0.0022	0.1089 ± 0.0025	8.5 σ	1.13 / 1.40
7°/22°, EquipartitionLsOr	0.0395	0.0787 ± 0.0022	0.1059 ± 0.0025	8.2 σ	1.12 / 1.39
30°/17°, NonLocal	0.0452	0.0565 ± 0.0019	0.0714 ± 0.0021	5.2 σ	1.00 / 1.20
30°/17°, EquipartitionLsOr	0.0484	0.0523 ± 0.0018	0.0683 ± 0.0021	5.8 σ	0.99 / 1.18

With 4000 ions under EquipartitionLsOr, As 30 keV gives a tail of 0.057 static and 0.069 at 300 K at 7°/22° (amorphous 0.009). At 30°/17° it gives 0.020 static and 0.028 at 300 K (amorphous 0.015).

Over seeds 1 to 5 (4000 ions, NonLocal, 7°/22°), the tail is larger at 300 K than static for every seed and both ions:

• B: 0.075-0.085 static, 0.105-0.118 at 300 K.
• As: 0.057-0.060 static, 0.065-0.072 at 300 K.

The dRp ratios have a seed-to-seed standard deviation of:

• B: 0.018 static, 0.037 at 300 K.
• As: 0.052 static, 0.047 at 300 K.

For As the dRp ratio does not rise with vibration (1.88 static, 1.80 at 300 K). More ions reach the tail, but the deepest ones stop shallower: the 99th percentile is 112 nm static and 98 nm at 300 K.

The local Oen-Robinson half does not remove the rise, so the rise is robust to this change of loss treatment (#225, part 3). The cause is not established: EquipartitionLsOr still keeps half of the nonlocal loss, so both treatments share the amorphous-average nonlocal approximation, and this comparison cannot rule it out as a contributor. Correlated vibration is absent: neighbouring atoms are displaced independently (“Thermal vibration”). No source read for #225 gives the sign of the effect of correlations on this tail, so its expected direction is not stated.

Comparison with the literature (#225)

Searched: measured depth profiles of B or As implanted into crystalline (100) Si near 7° tilt at room temperature, and measurements against tilt, twist or wafer temperature.

Not opened. The leads from the issue could not be read:

• Miyake et al., J. Electrochem. Soc. 130, 716 (1983), doi:10.1149/1.2119789. Paywalled at the publisher.
• Tian et al., Mikrochim. Acta (1992), doi:10.1007/BF01244469. Closed.
• Klein et al., IEEE Trans. Electron Devices 39 (1992), doi:10.1109/16.141226. Closed. It is a UT-MARLOWE paper, Tier C.

Several further candidates were not readable either:

• Semicond. Sci. Technol. 5 (1990), doi:10.1088/0268-1242/5/10/001, “Boron implants in <100> silicon at tilt angles of 0 degrees and 7 degrees”. The publisher served a bot check.
• J. Electrochem. Soc. papers on B channeling against tilt angle: doi:10.1149/1.2108785 and doi:10.1149/1.2085934. The publisher served a bot check.
• K. Nordlund, F. Djurabekova and G. Hobler, Phys. Rev. B 94, 214109 (2016), doi:10.1103/PhysRevB.94.214109. The publisher refused access and the repository copy did not respond.

None of these is cited for content here.

Read: measured profiles, not usable at matched conditions. D. Cai, N. Grønbech-Jensen, C. M. Snell and K. M. Beardmore, Phys. Rev. B 54, 17147 (1996), arXiv:physics/9901056 (read in full). The paper shows SIMS profiles of B in (100) Si:

• at 15, 35 and 80 keV, at tilt 7° and rotation 30° (Fig. 3);
• at 5 keV, at tilt 7° and rotation 7° (Fig. 4; the SIMS data are from K. Gärtner, M. Nitschke and W. Eckstein, Nucl. Instrum. Methods B 83, 87 (1993), which was not found openly readable).

It also shows As at 8°/30° (Fig. 6). These profiles were not digitized, for three reasons:

1. The paper does not give the in-plane reference direction of “rotation”, so the orientation cannot be mapped through docs/crystal-orientation.md.
2. It does not give the beam divergence, and it does not give the oxide or dose of the 5 keV measurement. Its 16 Å native oxide is a simulation setting.
3. Each measured line overlaps the UT-MARLOWE curves (Tier C) in the same vector figure.

A single room-temperature profile would also only show the model at 300 K. It could not tell the static and the vibrating lattice apart.

Read: a published simulation of the same question. M. I. Bratchenko, A. S. Bakai and S. V. Dyuldya, J. Phys. Stud. 13, 1601 (2009), doi:10.30970/jps.13.1601 (read in full). This is molecular dynamics with the authors’ own code, not a measurement. It simulates 15 keV B and As into (001) Si at 7° tilt and 30° rotation, for a static lattice and at 300 K, with uncorrelated Debye displacements. The paper reports:

• In the static lattice, the deeply channeled states are almost empty. The authors attribute this to blocking, which forbids capture into them by a single strong collision.
• At 300 K “thermal vibrations facilitate the volume capture of ions into the stable channeling mode”, that is, feeding-in.
• For As this lengthens the channeling tail substantially compared with the static lattice. For B the temperature effect on the tail length is “much weaker”.
• “As long as a target is ordered there is no chance to obtain truly random-equivalent doping profiles”.

Its introduction also states, citing SIMS work not read here, that off-axis SIMS profiles show long tails. It says the fraction of ions beyond the amorphous Gaussian reaches about 20-30 % below 20 keV.

The sign of the population effect here agrees with that paper: vibration feeds ions into channels at an off-axis direction. But the energy (5 and 30 keV here, against 15 keV) and the orientation (7°/22° here, against 7°/30°) differ, and so does the measure. The paper discusses the tail length, and for As at 30 keV our deepest percentile shortens. The paper uses the same uncorrelated vibration model, so it does not test that assumption. The off-axis rise is therefore consistent with one published simulation and not validated against measurement. Validating it needs a measured profile at matched, fully stated conditions, ideally at two wafer temperatures (docs/validation.md, “Channeling (M2)”).

The 7°/22° criterion (#225)

#180 asked that dRp at 7°/22° be within 15 % of amorphous. That bound is withdrawn. The beam lies 2.6° from a {100} plane and 2.7° from a {110} plane (7.0° from <100>), so it is a near-planar direction, not a random one. A tail of a few per cent of channeled ions is expected there, static or vibrating, and the second moment weights it heavily. Without a matched measurement, the 7°/22° runs are recorded-value regression checks (near_planar_7_22_static_lattice in tests/crystal_bca.rs and near_planar_7_22_at_300_k in tests/crystal_thermal.rs):

• The dRp ratio stays within 0.12 (B) or 0.16 (As) of the table value. That is three seed-to-seed standard deviations, rounded up.
• The 90th percentile stays within 30 % of the amorphous one.
• Static lattice: Rp stays within 10 %.
• At 300 K: the tail exceeds the amorphous tail by more than five binomial standard errors.

These are not agreement with experiment. They catch unintended changes of the model.

At 30°/17° the static lattice keeps both random-direction bounds. At 300 K the dRp ratio, 1.18 (B) and 1.17 (As), is above the 15 % bound. This is accepted as a property of the vibrating model, for three reasons:

• The bulk of the profile still matches: Rp and the 90th percentile are within 10 % of amorphous.
• The excess is the thermally fed tail measured above, at 5 σ.
• The one published study read finds that thermal vibration causes this kind of feeding-in.

The 300 K check asserts Rp and the 90th percentile within 10 %, and dRp below its recorded value plus 0.12 (random_direction_rp_at_300_k_matches_amorphous).

Unless marked EquipartitionLsOr, the ranges above use the nonlocal Lindhard-Scharff loss only (ElectronicLoss::NonLocal, the default), so a channeled ion loses as much electronic energy per unit path as in the amorphous target. With the local Oen-Robinson half (next section) a channeled ion loses less per unit path than one in a random direction. The off-axis tails under that mode are in the 16000-ion table above.

Local electronic loss (Oen-Robinson)

Under ElectronicLoss::EquipartitionLsOr half of the Lindhard-Scharff loss is nonlocal, taken along every flight segment (above). The other half is local: every partner of a collision step, so every lattice site that the search accepts within p_max of the path, takes the Oen-Robinson loss (OenRobinson::local_loss) at its own distance of closest approach r_min. That distance comes from the partner’s impact parameter, which with thermal vibration is the distance to the displaced position. The energy is the energy before the collision step. The losses Q_i enter the energy-conserving scaling above. O. S. Oen and M. T. Robinson, Nucl. Instrum. Methods 132, 647 (1976).

What is not modelled:

• The p_max cutoff on the local half. Sites farther than p_max from the path are not partners, so they take no local loss. The local half is therefore below half of Lindhard-Scharff wherever p_max is not large compared with the decay length a / 0.3 of the Oen-Robinson loss. The amorphous model has the same cutoff. In a channel, where the ion stays farther than p_max from most rows, this cutoff is what lowers the loss.
• Vibration beyond the displaced positions. The local loss depends on the vibration amplitude only through the displaced site positions; there is no other thermal term.
• Verified constants. The Oen-Robinson constants are not verified against the paper (OR_CONSTANTS_UNVERIFIED). CrystalMetadata::electronic_constants_unverified records this in the run metadata.

Measured by tests/crystal_electronic.rs (#[ignore]d; release build): R = E_local / E_nonlocal of the primary (recoils not followed), counted until its energy first falls below 80 % of the beam energy. Ar 20 keV into Si, static lattice unless noted, seed 1, the crystal’s p_max set to the amorphous constant-path radius 1.532 Å. Errors are delta-method standard errors.

Case	Histories	R	Notes
amorphous	4000	0.8645 ± 0.0026	the same at 30°/17° and 7°/22°
30°/17°	4000	0.9304 ± 0.0025	1.076 ± 0.004 × amorphous
30°/17°, 300 K	4000	0.9270 ± 0.0025	1.072 ± 0.004 × amorphous
7°/22°	4000	0.8744 ± 0.0030	1.011 ± 0.005 × amorphous
30°/17°, p_max = nn	4000	1.0033 ± 0.0026
along <110>	8000	0.2448 ± 0.0022	0.263 × the 30°/17° crystal value (0.9316, same n)
along <100>	8000	0.5600 ± 0.0010	0.601 × the 30°/17° crystal value
along <110> / <100>, 300 K	8000	0.2548 / 0.5719	30°/17° at 300 K: 0.9283

Along both axes the local share is far below the off-axis direction (by more than 150 standard errors), as channeling requires.

One direction is not a random medium (#250)

Along the fixed direction 30°/17° the crystal ratio is 7.6 % above the amorphous one at the same p_max (5.5 % at 63°/8°). #226 asked for 5 %. This is not a counting error. It is what the rule of angular averages of J. Lindhard, Mat. Fys. Medd. Dan. Vid. Selsk. 34, no. 14 (1965), section 5, eqs. (5.1)-(5.3), predicts for one direction. The rule states that a quantity linear in the probability of a particle being at a point of the crystal “has the same angular average as in a random system”, a system “with the same density, but without directional effects” (the paper names electronic and nuclear stopping at a fixed energy as examples, and calls this “compensation” of the directional effect). Only the average over directions equals the random medium. The channeling directions lie far below it (the <110> and <100> rows above), so the directions outside the channels lie above it.

The engine obeys the rule. direction_averaged_ratio_matches_amorphous averages over 1181 beam directions that cover, evenly in solid angle, the sixth of the sphere nearest to [100] (16 histories each, 18896 in all, a different seed per direction). It counts over a fixed path of 100 Å, because the paper limits the rule to small energy loss: the energy window above lets a channeled ion travel farther, so it weights the channels more.

Ar 20 keV into Si, 100 Å path window	R	× amorphous
amorphous, all directions	0.8416 ± 0.0011
crystal, all directions	0.8333 ± 0.0016	0.990 ± 0.002
crystal, all directions, 300 K	0.8365 ± 0.0016	0.994 ± 0.002
amorphous, 30°/17° (4000 histories)	0.8383 ± 0.0024
crystal, 30°/17° (4000 histories)	0.9120 ± 0.0021	1.088 ± 0.004

The direction average meets the 5 % of #226 (the test asserts max(3 σ, 5 %)). Its remaining 1 % is not asserted to be zero: the rule holds for energy-conserving motion, and the ions here slow down. Single directions of that run (16 histories, so noisy) give R from 0.37, near a <110> axis, to 1.09. Their median is 0.87, above the amorphous value.

fixed_direction_partners_lean_to_small_impact_parameters shows where the excess along 30°/17° comes from (BcaTally::partner reports every partner’s impact parameter p):

• The crystal makes as many collisions per unit path as the amorphous model (1.005 times).
• A partner at a given p takes the same local loss in both, within 3 % in each tenth of the disc area. The two engines share the loss function.
• The crystal’s partners are not uniform over the disc of radius p_max. The inner three tenths of its area hold 1.118 ± 0.004 of their uniform share and the outer three tenths 0.903 ± 0.004. The amorphous disc is uniform, and so is the crystal’s in the direction average (every tenth within 1.7 % of its share).

So along this direction the ions pass closer to the atoms than a uniform flux would. That is the excess of R, because the Oen-Robinson loss falls steeply with p.

Each mechanism of the collision step that could bias the partners was switched off in turn. The switches were temporary and are not in the tree. Same run as the table, 30°/17°, static:

• No deflection (the direction never changes; 100 Å path window): the disc is uniform and R is 0.8476, 1.011 × amorphous. So the lattice search and the count per path are right, and the excess needs the deflections.
• No exclusion of the sites just hit: R is 0.9297 instead of 0.9304.
• One partner per step only: R is unchanged (0.9304). Only 0.06 % of the steps have two partners at this p_max.
• The same sites at randomised positions (each site moved by its own fixed vector, uniform over a cube of edge a, which removes the short- range order and keeps the algorithm): R is the same with and without deflection (0.8645 and 0.8633, 100 Å path window), and every tenth of the disc is within 2.5 % of its share. So the half-tube search that starts at the collision point and the ordering by path distance do not bias the partners. The half-tube start was not switched off by itself. This medium gave 2 % more partners per path than the amorphous model, with and without deflection; that was not followed up.
• p_max: the excess of the lattice over the randomised sites hardly depends on it: 1.056, 1.053 and 1.046 at p_max = 1.53, 2.35 and 3.2 Å (energy window). It is not an effect of the cutoff.

No repeated collision with one site was found in the 30°/17° run (0 of 137251 partners). Thermal vibration up to 3000 K and q_max hardly change the excess (#250).

The energy window is not neutral either. In the run without deflection it gives R = 0.9127 and uneven tenths, where the path window gives 0.8476 and a uniform disc. A window that ends a history on a large loss selects histories when successive collisions are correlated, as they are along a fixed line through a lattice. With deflection the two windows give 1.076 and 1.088 × amorphous.

The fixed-direction quotients are therefore recorded values (fixed_direction_local_to_nonlocal_ratio_is_above_amorphous: 1.076 static and 1.072 at 300 K, each within 0.02), not agreement with a reference. A run along one off-axis direction takes a few per cent more local electronic loss than the amorphous model at the same p_max.

Thermal vibration (step 21b)

With CrystalTarget::thermal set, each lattice site the particle meets is displaced from its static position by a random vector u whose three Cartesian components are independent Gaussians of standard deviation u1, the one-dimensional RMS amplitude of the Debye model (crate::ion::crystal::debye) for the site’s species at the target temperature Thermal::temperature_k. Displacements of different sites are uncorrelated; a correlated mode is not implemented. The issue that asked for this model (#181) names D. S. Gemmell, Rev. Mod. Phys. 46, 129 (1974), doi:10.1103/RevModPhys.46.129, as its reference; that paper could not be opened, so nothing here is taken from it and none of its equations is cited. The amplitude and its source are those of the debye module.

• Encounter. A site’s displacement is drawn when the site first appears in a search of the current flight line (a straight path of one particle in one crystal) and is kept while the line goes on: the next search segments of a channel see the same displaced atom. Any collision ends the line (and so does a change of direction at a surface or leaving the crystal), so every site met afterwards, by the projectile or by a recoil, gets a fresh draw. A site is therefore displaced once per encounter. Correlations in time (an atom met again later by the same cascade) are not kept.
• What the displacement changes. The impact parameter, the path distance (so the order of the partners and the simultaneity test), the direction n_i from the line to the atom, and so the deflection, the recoil direction and (under EquipartitionLsOr) the local electronic loss. A site’s region (and so whether it is a partner at all) and the starting point of its recoil are its static lattice position, so that recoils, vacancies and tallies stay on lattice sites and a recoil’s own site is recognised as such; the difference is u1, 0.065 Å for Si at 300 K (THETA_D_SI), against 2.35 Å between neighbours.
• Search margin. The static lattice is searched with the radius p_max + m over the path range [-m, L + q_max + m], with m = THERMAL_MARGIN_U1 * u1 (largest u1 of the crystal), and the static selection rules (0 <= s <= L + q_max, p <= p_max) are then applied to the displaced positions. A site whose displacement perpendicular to the path exceeds m can be missed: the probability is exp(-m^2 / (2 u1^2)) = exp(-18), about 1.5e-8, per site.
• Random numbers. The displacements come from the history’s counter-based stream (crate::rng::stream), from a copy positioned at word 2^66, a segment that neither the transport draws nor the lattice translation (word 2^67) reach. Each displacement consumes exactly six words (sample_gaussian_3d), drawn in the documented order of the search results (crate::ion::crystal::search, “Ordering”), and the primary and its recoils share the stream in the order they are followed. So a history’s result depends only on (seed, index), never on the thread count, and turning vibration on changes no other draw.
• Static limit. With every amplitude zero (T = 0 and Thermal::include_zero_point = false, a test switch) the thermal code path gives the static model bit for bit (tests/crystal_thermal.rs). With the zero-point term (the physical case) T = 0 still vibrates: u1 = 0.045 Å for Si (0.088 Å at 600 K).
• Not changed. The lattice constant is the cited one at its own temperature (no thermal expansion). The electronic-loss model is the one of the static lattice: the nonlocal loss along a segment, and under EquipartitionLsOr the local loss at each partner’s closest approach. The local loss sees the vibration only through the displaced position (the impact parameter, so r_min); it has no other thermal term.

DISPLATH (Tier A, above; same commit) also perturbs lattice atoms by independent Gaussian displacements of a Debye amplitude, but draws them once per cell and cascade (LatticeSiteCoordinatesBase! and TemperatureToSigma in src/geometry.jl); the per-encounter scheme here is this crate’s, as #181 asks, and no DISPLATH code was ported for it.

What is not here

The lattice is perfect (static, or vibrating as above) and stays perfect: a displaced atom leaves no vacancy, so a later particle can still collide with its site (dynamic damage, dechanneling and amorphization are later steps). The impact-parameter-dependent electronic loss is the local Oen-Robinson half, taken only at sites within p_max (see “Local electronic loss” for what that leaves out). Boundaries of a periodic voxel axis do not continue the lattice. The equation-level treatment of Robinson and Torrens, Phys. Rev. B 9, 5008 (1974) (MARLOWE, Tier C) was not available to read when this was written; the choices above are DISPLATH’s, not quotations of that paper.

Structs

CrystalMetadata

Run metadata of one crystal of a Bca (Bca::crystal_metadata): what a result needs to be reproduced from its own header.

CrystalTarget

A crystal for Bca::with_crystal: the lattice, the orientation of the target in the lab frame and the three search parameters.

Thermal

Thermal vibration of the lattice atoms of a CrystalTarget (step 21b; module docs, “Thermal vibration”).

ThermalMetadata

The thermal part of CrystalMetadata.

Constants

DEFAULT_Q_MAX_FRACTION

Default q_max, as a fraction of the nearest-neighbour distance (module docs, “Choice of the search parameters”).

THERMAL_MARGIN_U1

Search margin of the thermal model, in units of the largest one-dimensional RMS displacement u1 of the crystal’s species (module docs, “Thermal vibration”).

THERMAL_SAMPLING

How Thermal displacements are sampled, as recorded in ThermalMetadata::sampling.