## 2 Experimental Techniques and Methodology
For each pattern, data were recorded in the angular range \(5^{\circ}\leq 20\leq 121^{\circ}\) with a step size of 0.05\({}^{\circ}\) using neutrons of wavelength 1.2253 A. The X-ray measurements were done on a Bruker D8 Advance apparatus working in the Bragg-Brentano geometry in the \(\theta\)-\(\theta\) mode with Cu K\(\alpha\) radiation and a graphite monochromator in the diffracted beam. The usual conditions for the measurements are 5-120\({}^{\circ}\) (20), step 0.04\({}^{\circ}\), time per step 20-60 s.  
To determine site occupancies, the combined Rietveld analysis of XRD and ND data was conducted using the FULLPROF program in the multipattern mode (joint refinement). For refinement, the background was either interpolated between the peaks or refined with a polynomial. The pseudo-Voigt function was used to define the peak shape in both XRD and ND patterns. Because Re and Co present a  
Figure 3: Diagrams showing as parallel lines the compositions of a site occupied by the three atoms Mo, Re, and Ni (**A**) or Cr (**B**), giving the same diffraction scattering density with X-rays (in blue; it is assumed that the scattering is proportional to the average \(Z\)) and neutrons (red). non-negligible absorption coefficient for neutrons, an absorption correction was also applied by calculating the values of the absorption coefficient (\(\mu r\)) for each composition. To calculate the absorption coefficient, the density of non-compacted powders was estimated to be \(\sim\)40% of the density of the material. Also, micro-absorption for the flat-plate geometry of X-ray measurements was taken into account in particular in order to be able to describe both XRD and ND data with a unique structural model including unique displacement parameters. The site occupancies are refined using a composition constraint and maintaining full occupancy at each site. The way to handle such complex relations between the parameters is described in [5]. First-principles calculation of stable and metastable intermetallic phases has become routine. However, systematic calculations in ternary systems may become rather complex. For a crystal structure with \(n\) independent sites in the unit cell, the ordered distribution of three atoms among each site generates \(3^{n}\) configurations. If the calculation is feasible for intermetallic compounds with two or three sites like C15 and C14, it gets much more complex for phases like \(\chi\) and \(\sigma\) with four and five sites generating \(3^{4}=81\) and \(3^{5}=243\) ordered configurations. For this reason, the ZenGen code was developed, allowing us to perform systematic calculations of all the configurations generated by the distribution of the three elements among the sublattices [6]. Above six sites, the computation time may become prohibitive due to the increase in both the number of configurations and the complexity of the structure; simplifications should be adopted, such as merging different crystallographic sites. The presented calculations are based on the density functional theory (DFT) and were carried out using the VASP package [7, 8]. The generalized gradient approximation (GGA) was used with the Perdew-Burke-Ernzerhof (PBE) exchange and correlation energy functional [9, 10]. An energy cutoff of 400 eV was used for the plane wave basis set, and a dense grid of \(k\)-points was used with the sampling generated by the Monkhorst-Pack procedure [11]. For each ordered configuration, both the internal atomic coordinates and the lattice parameters were fully relaxed using successive steps, keeping the original cell structure. The conditions for the calculations (method, accuracy) have been detailed previously [3, 4]. All the calculations have been handled using the ZenGen code [6] (generation of input files, job monitoring, extraction of results). This procedure allows the calculation of the heat of formation of every compound generated by assigning each element to each crystal site, and gives access to the stability of the phase.