To model the experimental system, we constructed slabs of 4\(\times \) 4\(\times \) 2 size formed by four layers (192 atoms each). The relaxed atomic positions are represented in Fig. REF , together with a lateral view of the spin isosurfaces discussed in the main text. All the structures were fully optimized without constraints until the forces on each atom were smaller than 10\(^3\)  Ry/au and the energy difference between two consecutive relaxation steps less than 10\(^4\)  Ry. The Brillouin zone was sampled by a \(\Gamma \)  centered 3 \(\times \)  3 \(\times \)  1 k-point Monkhorst-Pack [1]} mesh for structural optimization and 6 \(\times \)  6 \(\times \)  2 for the self-consistent field (SCF) calculations. We built two supercells. The supercell A is formed by a single Fe impurity substituting a Nb atom and 63 Nb and 128 Se atoms. The supercell B includes a 10% random substitution of Se by S atoms, resulting in 1 Fe, 63 Nb, 115 Se and 13 S atoms.
