Scattering modes

IRMA describes a material's thermal response with three choices: which coherent-elastic treatment to use (iel), which elastic format to write (elastic_mode), and which inelastic engine to run (inelastic_mode). S(α,β) is the thermal scattering law in the dimensionless momentum and energy transfer; the theory page defines it and the conventions used here. The selectors live on different cards (the numbered records of the LEAPR-style input file; see the input file reference). iel is set on Card 5. When iel=10, Card 6b adds the other two selectors, elastic_mode and inelastic_mode, which can be combined in any way. This page describes each selector, shows how to read them off an input file, and gives guidance for picking the right combination for cubic versus noncubic and isotropic versus anisotropic materials.

The three axes at a glance

Axis	Card / field	What it controls	Values
iel	Card 5, field 4	Coherent-elastic (Bragg) treatment	0 none/incoherent; 1–6 built-in materials; 10 generalized (any crystal)
elastic_mode	Card 6b, field 1 (iel=10 only)	ENDF elastic format written	1 = SEF (single-channel elastic format); 2 = MEF
inelastic_mode	Card 6b, field 4 (iel=10 only)	Inelastic S(α,β) engine	0 phonon expansion from a DOS; 1 directional incoherent; 2 coherent one-phonon

elastic_mode and inelastic_mode live on Card 6b, which is present only when iel=10. For the built-in materials (iel=1–6) and for iel=0, IRMA uses the classic kernels and the two Card 6b selectors do not apply.

Axis 1 (iel): coherent-elastic treatment

iel (Card 5, field 4) selects how Bragg (coherent-elastic) scattering is handled.

iel	Treatment	Notes
0	None / incoherent	No coherent elastic; with no translational mode IRMA flags incoherent elastic
1–6	Built-in coherent elastic	Graphite, Be, BeO, Al, Pb, Fe (in that order)
10	Generalized	Bragg edges computed from the crystal structure you supply on Cards 6c–6d; works for any material

The built-in options reproduce the crystal tables of LEAPR's historical coher routine. iel=10 is recommended for new evaluations: it computes Bragg edges directly from your unit cell and unlocks the elastic_mode and inelastic_mode axes.

The generalized path also changes what the Card 4 za must contain. With iel=10, za must encode the physical nuclide as 1000·Z + A (with A = 0 for the natural element) so that the principal scatterer can be matched to a Card 6d atom type. Some thermal libraries use a TSL-convention ZA instead; that convention will not match here.

Axis 2 (elastic_mode): which elastic format to write

elastic_mode (Card 6b, field 1) is only meaningful when iel=10. It selects the ENDF elastic representation and therefore the LTHR flag your downstream codes will see.

elastic_mode	Name	What is written	Resulting LTHR
1	SEF (single-channel elastic format)	One elastic component per tape (each species' tape), scaled to carry the total bound elastic scattering (mechanism below)	1 (coherent) or 2 (incoherent), per species
2	MEF (Mixed Elastic Format)	Every atom gets both a coherent (per-atom Bragg edges) and an incoherent part	3 (mixed)

For a single-atom material under SEF, the component with the larger bound cross section, coherent or incoherent, is kept and scaled to carry the total bound scattering cross section. For a polyatomic cell, SEF selects one designated-coherent (DC) atom, the one minimizing the incoherent contribution; that species' tape (ENDF output file) carries the coherent elastic, and the others carry incoherent elastic with a redistribution factor. The scaling has a physical consequence for coherent-kept materials: coherent elastic scattering does not exist below the first Bragg edge, so the SEF elastic cross section is zero there. MEF carries the incoherent component separately and does not have this gap.

MEF needs downstream support: elastic_mode=2 writes LTHR=3, so confirm that your transport and processing chain understands the mixed elastic format before choosing it. When in doubt, SEF (elastic_mode=1) is the conservative, widely supported choice.

Axis 3 (inelastic_mode): the S(α,β) engine

inelastic_mode (Card 6b, field 4, iel=10 only) chooses how the inelastic S(α,β) is built.

inelastic_mode	Engine	Debye-Waller	Needs phonopy?	Reads the DOS cards?
0	Phonon expansion from a tabulated DOS	Isotropic (scalar)	No	Yes (Cards 11–19, with Card 12e)
1	Directional incoherent approximation	Directional	Yes	No
2	Exact coherent + incoherent one-phonon, incoherent multiphonon	Directional	Yes	No

inelastic_mode=0 is the classic LEAPR physics level: an isotropic Debye-Waller factor and a phonon expansion built from a scalar phonon density of states given in the input file, with no phonopy required. inelastic_mode=1 computes a directional Debye-Waller factor for the coherent elastic and an S(α,β) computed in-process from the phonopy calculation, using an incoherent-approximation one-phonon term plus incoherent-approximation multiphonons. inelastic_mode=2 is the same as mode 1 except that the one-phonon term is exact (coherent plus incoherent) on top of the incoherent-approximation multiphonons.

In modes 1 and 2 the one-phonon term comes from the phonon modes of the phonopy q-mesh, so it is zero for energy transfers below the lowest kept mesh-mode energy; a finer mesh reaches lower. The run prints that energy and the number of output energies below it.

Even in mode 2, only the one-phonon term is coherent. The multiphonon orders (n ≥ 2) use the incoherent-approximation model (a per-atom self kernel scaled by the atom's total scattering cross section), so coherent interference is dropped in the tail. The engine records multiphonon_model = "incoherent_approximation" in its run metadata (the standalone engine's JSON output and the spectra results' engine metadata); ENDF tapes do not carry this label.

Modes 1 and 2 ignore the DOS cards. inelastic_mode=1/2 build the inelastic section (MF7/MT4) and the Debye-Waller factors entirely from the phonopy calculation: the continuous-DOS, translational, and oscillator detail cards (Cards 11–19) are not read, and Card 12e partial spectra are rejected (nspec must be 0). The input file supplies only the temperature cards; everything else comes from phonopy. The mixed-moderator (nss>0), cold-hydrogen (ncold), and Sköld (nsk) options are also unavailable with modes 1/2 and are rejected at parse time.

What modes 1 and 2 require

When inelastic_mode=1 or 2, supply the phonopy mesh and sampling controls between Cards 6d and 7:

• Card 6f: the phonopy.yaml path, then mesh_nx mesh_ny mesh_nz ncpu use_born (and a BORN path if use_born=1). Force constants are read from the yaml if embedded, otherwise discovered next to it (force_constants.hdf5, FORCE_CONSTANTS, FORCE_SETS, in that order). The process working directory is never consulted.
• Card 6g: the powder-averaging sampling controls (next section).

These modes need the optional phonopy dependency:

pip install -e ".[phonopy]"   # + the phonopy-backed modes (inelastic_mode=1/2)

Choosing a combination

Cubic vs noncubic materials

For a cubic material the Debye-Waller factor is isotropic by symmetry, so inelastic_mode=0 (with a good DOS) captures the physics and is fast. For a noncubic (anisotropic) crystal, the directional engines (1/2) are the point of iel=10: they carry the full Debye-Waller tensor rather than collapsing it to a scalar.

Isotropic vs directional Debye-Waller: why it matters

The Debye-Waller factor suppresses one-phonon scattering exponentially with momentum transfer Q. In an anisotropic crystal that suppression is direction-dependent, and the powder average of the exact directional factor is not the same as applying a single orientation-averaged (isotropic) factor.

In graphite, for example, the perpendicular and in-plane Debye-Waller terms differ by roughly a factor of 6.6 (W_c / W_ab ≈ 6.6). An isotropic treatment exponentially over-suppresses the one-phonon S(α,β) at high Q relative to the exact directional powder average. IRMA's inelastic_mode=2 performs the directional powder average exactly; inelastic_mode=0 (isotropic Debye-Waller) on the same phonon calculation shows the over-suppressed roll-off. Mode 0 also differs from mode 2 in the one-phonon term itself, which it builds in the incoherent approximation from a scalar DOS.

The exact one-phonon treatment is cross-validated against Euphonic, an independent code run on the same phonon calculation, for graphite, Be, and BeO: coherent component against coherent component, the integrals of the symmetric S(α,β) agree to ratios of 1.00001 (graphite), 1.0002 (Be), and 0.9998 (BeO). For strongly anisotropic crystals, prefer inelastic_mode=2; reserve inelastic_mode=0 for cubic materials or for reproducing existing isotropic evaluations.

Card 6g: sampling controls

Card 6g (modes 1/2 only) sets the powder-averaging resolution and the multiphonon order:

ndir  mpdir  [auto_order]  /

Field	Name	Meaning
1	ndir	One-phonon powder-average directions, for the coherent and the incoherent term (golden-spiral quadrature). Must be ≥ 1
2	mpdir	Multiphonon powder-average directions. Converges by about 50–100 directions; 1000 is recommended, and the cost grows linearly with mpdir. Must be ≥ 1
3	auto_order (optional)	0 = honor the Card 3 nphon multiphonon order verbatim (default); 1 = auto-size the multiphonon order

The incoherent powder average is always the exact numerical orientational average; there is no method selector.

Recommended Card 6g

The recommended production setting (the GUI and irma mlip emit default) is:

10000 1000 1 /

That is 10000 one-phonon powder directions, 1000 multiphonon directions, and auto-sized multiphonon order: the same sampling the validation campaign ran. The powder averages converge well below these counts, so a faster exploratory run at 4000 200 1 is fine; the default carries the cross-code-comparison headroom.

Auto-sizing is what keeps the tabulated S(α,β) converged at high Q. With auto_order=1, IRMA raises the multiphonon order as needed to converge the anisotropic-Debye-Waller Poisson sum at the grid's Q_max: graphite on a grid reaching Q_max ≈ 100 1/Å needs order ≈ 223, so a fixed nphon=100 truncates the table at high Q, and the engine warns when this happens. Leaving auto_order=0 honors the Card 3 nphon exactly. That is appropriate when you want full control, but make sure nphon (or auto_order=1) gives the tabulated S(α,β) genuine support over your requested beta grid. Energy transfers beyond the tabulated range are covered downstream by the short-collision-time extension of THERMR, NJOY's thermal processing module.

Bragg-edge grouping (optional)

For materials with many closely spaced high-energy Bragg edges, Card 6b accepts two optional trailing fields that enable ENDF-102 §7.2.2 Bragg-edge grouping:

elastic_mode  nat  nspec  inelastic_mode  [bins_per_decade]  [threshold_eV]  /

Field	Name	Meaning
5	bins_per_decade	Number of log-uniform bins per decade above the threshold. 0 or absent = grouping off (default)
6	threshold_eV	Energy above which edges are merged. Default 1 eV

Above threshold_eV, the dense edge steps are merged into bins_per_decade log-uniform bins per decade with structure-factor-weighted placement. Grouping keeps the cumulative S exact at every bin boundary and every temperature. Within a bin the edges become one step, so the pointwise cross section changes slightly, and more at temperatures away from the first; the run prints the per-temperature error. It is off by default.

Grouping merges only the edges above the threshold and keeps every sub-threshold edge raw, whereas the ungrouped writer applies NJOY's tolerance thinning to the whole list. A very low threshold_eV on a material whose Debye-Waller factor suppresses the high-energy edges can therefore produce a larger tape than grouping off. Use the default 1 eV threshold unless the dense structure you want compressed actually sits below it.

Putting it together: a mode-2 input file

A minimal iel=10 / inelastic_mode=2 Card 5–6g block looks like this (graphite, exact one-phonon, auto-sized multiphonon order):

 11.898 4.739 1 10 0 0 /       --- Card 5: awr, spr (free-atom sigma), npr, iel=10, ncold=0, nsk=0
 0 /                            --- Card 6: nss=0 (no secondary scatterer)
 1 1 0 2 /                      --- Card 6b: SEF, nat=1, nspec=0, inelastic_mode=2
 2.464 2.464 6.711 90 90 120 / --- Card 6c: hexagonal graphite lattice
 6 12 11.898 6.6484 0.0 4 /     --- Card 6d: C-12, b_coh [fm], sigma_inc [b], 4 positions
 0 0 0.25  0 0 0.75  ...  /     --- Card 6d: fractional positions
 'phonopy.yaml' /               --- Card 6f-1: phonopy.yaml path
 40 40 40 14 0 /                --- Card 6f-2: mesh, ncpu=14, use_born=0
 10000 1000 1 /                 --- Card 6g: ndir, mpdir, auto_order=1

The numbers in this snippet are illustrative: the lattice constants, cross sections, and positions show the input file shape, not a validated evaluation. The production-ready input file is examples/tsl/graphite_mode2.input, a complete, runnable mode-2 input file at the recommended production sampling with auto-sized multiphonon order (see examples/tsl/README.md for the run command). The tests/mode2_euphonic_n1_validation/*.template files are the cross-code validation harness: they deliberately compute the one-phonon term only and are not production evaluations.