Disorder-induced damping of spin excitations in Cr-doped BaFe$_2$As$_2$

TL;DR

This study demonstrates that disorder from Cr doping in BaFe2As2 significantly dampens and localizes spin excitations, overshadowing charge doping effects, as shown by RIXS and ARPES measurements and supported by a localized spin model.

Contribution

It introduces a disorder-inclusive localized spin model that explains the damping and localization of spin excitations in Cr-doped BaFe2As2, highlighting disorder's dominant role.

Findings

Fe-derived magnetic excitations are overdamped at high Cr doping

Cr doping increases electronic localization and suppresses nematic band splitting

Disorder effects dominate over charge doping in this Hund's metal

Abstract

In doped Hund's metals, such as the iron-based superconductors, effects like charge doping and chemical pressure are often considered the dominant factors. Partial chemical substitution, however, inevitably introduces disorder. Here, we investigate spin excitations in Ba(Fe$_{1-x}$Cr$_x$)$_2$As$_{2}$ (CrBFA) by high-resolution resonant inelastic x-ray scattering (RIXS) for samples with $x = 0, 0.035,$ and $ 0.085$. In CrBFA, Cr acts as a hole dopant, but also introduces localized spins that compete with Fe-derived magnetic excitations. We found that the Fe-derived magnetic excitations are softened and damped, becoming overdamped for $x = 0.085$. At this doping level, complementary angle-resolved photoemission spectroscopy measurements (ARPES) show increased electronic localization and a suppression of the nematic $d_{xz}/d_{yz}$ band splitting present in the parent compound. We thus propose a localized spin model that explicitly incorporates substitutional disorder and Cr local moments, successfully reproducing our key observations. Our findings reveal a case where disorder dominates over charge doping in the case of a Hund's metal.

Figures

Tables

Table 1. Table S1: Coupling constants for exchange interactions between next-nearest neighbor sites in a square lattice, the ratio J 2 / J 1 J_{2}/J_{1} depends on the type of bond formed by Fe-sites and/or Cr-sites.

$⟨⟨i,j⟩⟩$ bond type	$J_2/J_1$
Fe - Fe	2
Cr - Cr	0
Fe - Cr	1

Equations

$${\chi }^{\prime \prime }(\omega )=\frac{1}{1-e^{-\hslash \omega /k_{B}T}}\cdot \frac{I_0\gamma \omega }{({\omega }^2-{\omega }_0^2{)}^2+4(\gamma \omega {)}^2}$$

$$S_{\text{C}}^{\mu \nu }(\mathbf{q},\omega )=\frac{1}{2\pi N}\sum _{ij}{\int }_{-\infty }^{\infty }dt⟨S_i^{\mu }(t)S_j^{\nu }(0)⟩e^{\text{i}\omega t-\text{i}\mathbf{q}\cdot {\mathbf{R}}_{ij}}.$$

$${\bar{S}}_{\text{C}}^{\mu \nu }(\mathbf{q},\omega )=\frac{t_{\text{max}}}{2\pi }{\left[{⟨S_{\mathbf{q}}^{\mu }(\omega )S_{-\mathbf{q}}^{\nu }(-\omega )⟩}_{{\text{C}}_S(0)}\right]}_{D_x},$$

$$S_{\mathbf{q}}^{\mu }(\omega )=\frac{1}{\sqrt{N^{}{N}_t}}\sum _{i,n}{S}_i^{\mu }(t_n)\exp (\text{i}\omega t_n-\text{i}\mathbf{q}\cdot {\mathbf{r}}_i).$$

$$\frac{d{\vec{S}}_i(t)}{dt}={\vec{h}}_i(t)\times {\vec{S}}_i(t).$$

$${\mathcal{S}}^{\mu \nu }(\mathbf{q},\omega )=\frac{\beta \omega }{1-\exp (-\beta \omega )}{\bar{S}}_{\text{C}}^{\mu \nu }(\mathbf{q},\omega ).$$

$$I_{fluo}=\alpha \exp (-a\omega )\omega (1-g_{{\Gamma }_1})+\beta \exp (b\omega )g_{{\Gamma }_1}+\gamma \exp (c\omega )g_{{\Gamma }_2}$$

$$g_{{\Gamma }_{1,2}}=(\exp (-(\omega -{\omega }_{1.2})/{\Gamma }_{1,2})+1{)}^{-1}$$

Full text

Introduction: Electronic correlated materials often display multiple competing interactions, each associated with one or more energy scales. Partial chemical substitution tunes interaction energy scales to stabilize a desired ground state. This gives rise to rich composition ($x$) vs. temperature ($T$) phase diagrams, characteristic of correlated systems [1].

In iron-based superconductors (FeSCs) [2], chemical substitutions often induce the suppression of the antiferromagnetic (AFM) ordered phase and the emergence of high-temperature superconductivity (HTSC) [3, 4, 5, 6, 7, 8, 9]. While these effects are typically attributed to charge doping or chemical pressure [10], the role of disorder, an inevitable consequence of substitution, remains unclear. Early scenarios for the physics of FeSCs proposed that disorder could be the driving mechanism behind the observed $x$ vs. $T$ phase diagrams [11, 12] but subsequent approaches favored charge doping and chemical pressure effects [13, 14].

The parent compound BaFe2As2 exhibits a stripe-type AFM transition at $T_{\text{N}}\approx 134$ K, accompanied by a tetragonal-to-orthorhombic structural transition [15, 16]. Substitutions suppressing these transitions may or may not induce superconductivity (SC). For instance, Mn substitution (MnBFA) introduces minimal charge doping, mostly affecting the electron pockets [17, 18, 19], and some attention was devoted to the effect of disorder acting as the main tuning knob in this phase diagram [20, 21, 22, 23]. In contrast, Cr substitution (CrBFA) acts as a hole dopant [24], but the evolution of $T_{\text{N}}$ in CrBFA and MnBFA follows a similar trend as a function of $x$, with no emergence of SC for both substitutions. Therefore, charge doping alone cannot explain their phase diagrams.

Previous studies have shown that hole doping enhances electronic correlations in BaFe2As2 [25, 26], an effect often attributed to proximity to a Mott insulating state. However, experimental observations show that the ground state remains metallic even in regimes with localized magnetism. This apparent contradiction can be resolved by the concept of spin-orbital separation in multiorbital systems: the Hund’s coupling delays spin screening relative to orbital screening, leading to the so-called Hund’s metal behavior. BaFe2As2 is currently recognized as a Hund’s metal, where strong Hund’s coupling and orbital selectivity govern the degree of electronic correlations, which are expected to increase further with hole doping toward half-filling [27, 28, 29, 30]. While K substitution on the Ba site preserves itinerancy and superconductivity even at high hole doping, direct substitution at the Fe site, as in the case of Cr, introduces carriers that enhance localization and Fe-sublattice disorder. Resonant inelastic x-ray scattering (RIXS) is proven to be a powerful tool to investigate magnetism with the contrasting itinerant character of electrons and localized character of spins within FeSC phase diagrams [31, 29, 32, 33].

In this work, we investigate how spin excitations originating from Cr and Fe compete for the ground state properties across the $x$ versus $T$ phase diagram of CrBFA. To this end, we performed high-resolution RIXS, complemented by angle-resolved photoemission spectroscopy (ARPES) experiments and theoretical modeling. Our combined experimental and theoretical analysis shows that the key features of the data can be understood within a localized spin model, where disorder effects introduced by Cr moments play a central role. This reveals a case of a correlated metal in which disorder dominates over the charge doping effect.

Our results explain why the suppression of the AFM order in CrBFA and MnBFA does not follow the expected trend for charge doping, where Cr nominally introduces twice as many holes as Mn [34, 17]. They may also help explain the $x$–$T$ phase diagrams of Cr/Mn substituted CaKFe4As4, where the suppression of both HTSC and AFM order similarly deviates from charge doping expectations [35, 36]. Altogether, our findings suggest that the absence of SC in these systems is likely linked to the joint effect of increasing correlation and disorder effects on the Fe lattice introduced by Cr and Mn impurities.

Materials and methods: Ba(Fe1-xCrx)2As2 single crystals, with different concentrations of Cr, were grown using the In-flux method [37]. Sample characterization methods are described in the Supplemental Material S2 [38, 39, 40].

The $x=0,0.035,$ and $0.085$ samples (denominated as BFA, Cr$3.5\%$, and Cr$8.5\%$, respectively) were selected for RIXS experiments at the I21 beamline at Diamond Light Source [41], with a detected photon energy resolution of $\approx 30$ meV. The energy of the absorption spectra $L_{3}$ edge maximum at $709.5$ eV was chosen for the momentum-dependent studies. The methods of the sample preparation and data acquisition are detailed in the Supplemental Material S2.

We define the wavevector $\mathbf{q}$ in reciprocal lattice units (r.l.u.) as $(H,K)=(q_{x},q_{y})a/2\pi$, where $a=5.598$, 5.590, and 5.616 Å are the in-plane lattice parameters for BFA, Cr$3.5\%$, and Cr$8.5\%$, respectively. These correspond to the magnetic unit cell shown as a magenta dashed square in Fig. 1(a). Since the samples are twinned, we adopt $a=b$ and do not distinguish orthorhombic $a$ and $b$ directions. In this Brillouin Zone (BZ), we probe a maximum momentum of $|\mathbf{q}|\approx 0.56$ r.l.u., indicated by the blue circle in Fig. 1(b). High-symmetry points are labeled using the 1Fe BZ (black square) and consistently used for ARPES and theoretical results. ARPES acquisition details are as described in refs. [17, 24] and in the Supplemental Material S2.

To model the magnetic excitations, we employ a minimal antiferromagnetic (AFM) XXZ $J_{1}-J_{2}$ model on a square lattice, with anisotropy favoring the spins to lie on the $S_{x}-S_{y}$ plane. Classically, this system exhibits checkerboard (Néel) order for $J_{2}<J_{1}/2$ and stripe order for $J_{2}>J_{1}/2$. Cr doping ($x$) is treated as a local disorder in $J_{2}$, the coupling controlling the transition between the ordered phases. At $x=0$, the system has stripe order; at $x=1$, it displays Néel order. For intermediate doping, we randomly replace a fraction $x$ of Fe sites with Cr, assigning $J_{2}$ between sites $i$ and $j$ as follows: (i) Fe–Fe, $J_{2,ij}=J_{2}(x=0)$; (ii) Cr–Cr, $J_{2,ij}=J_{2}(x=1)$; (iii) Fe–Cr, $J_{2,ij}=[J_{2}(x=1)+J_{2}(x=0)]/2$. The classical ground state is determined via Monte Carlo simulations and annealing. Excitation spectra are computed by integrating the semiclassical equations of motion: $d\vec{S}_{i}/dt=\vec{h}_{i}\times\vec{S}_{i}$, where $h_{i}^{\alpha}=\sum_{j}J_{ij}^{\alpha}S_{j}^{\alpha}$, is the local exchange field at site $i$ and $\alpha=x,y,z$. Finally, we Fourier transform $\vec{S}_{i}(t)$ to obtain the dynamical spin structure factor for comparison with experiments. Details are provided in the Supplemental Material S1 [42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58].

Results and discussion: In Fig. 1(c)-(h), we show an overview of our RIXS momentum-dependent spectra for each sample. The probed reciprocal space high-symmetry directions and range are represented by the black dashed lines and blue circle in Fig. 1(b). A background subtraction was performed to better visualize and fit the quasi-elastic region (which also features low-energy phonon excitations) and magnon excitations. The subtraction method is explained in the Supplemental Material S2.1 [59, 60, 61].

The main magnon dispersion is clear even for the Cr$8.5\%$ substitution, for which the elastic line is enhanced and the magnon is softened. The magnon appears as a shoulder of the quasi-elastic peak at low $\mathbf{q}$ and is evident at maximum momentum transfer. To gain more quantitative insight, spectra were fitted to a model composed of an elastic (Gaussian), a phonon (Gaussian), a principal magnon peak, and a secondary inelastic magnonlike peak (Gaussian).

The main magnon peak can be fitted to the imaginary part of the dynamic susceptibility of a damped harmonic oscillator (DHO) $\chi^{\prime\prime}(\omega)$, as shown in Equation 1 [63, 64], where the $\mathbf{q}$ dependency reveals the magnon dispersion, $I_{0}$ is an intensity constant, $\gamma$ is the damping coefficient, and $\omega_{0}$ is the undamped (bare) frequency of the spin excitation as a function of $\mathbf{q}$. The Bose-Einstein distribution is included to reproduce thermal effects:

[TABLE]

We refer to the peak energy of the magnetic excitation as $\omega_{\text{peak}}$, i.e., the energy at $\chi^{\prime\prime}(\omega)$ reaches its maximum.

We note that it is necessary to include the second magnonlike peak in the fitting to accurately describe the shape and sharp width of the main magnon peak, especially for the highest $\mathbf{q}$ values and in the case of the parent compound. The inclusion of this peak is justified by the inspection of Fig. 2(a), wherein we show a direction and polarization survey of our results for the parent compound. As can be observed, a second excitation appears riding on top of the background fluorescence and peaking around $\approx 450-550$ meV. It is clearer for the $[100]$ $(\pi,\pi)$ direction, as indicated by the orange and blue arrows.

Figs. 2(b) and (c) show examples of fittings for the BFA and Cr$3.5\%$ samples. The data after the background subtraction are shown as blue dots, the total fitting as the red line, and the shaded regions as the different fitting contributions, including the elastic line, phonon, main magnon, and secondary magnon-like peak, shown as pink, green, blue, and orange areas, respectively. All the fittings for different momentum points for all samples and directions are shown in the Supplemental Material S2.1. The need for the secondary magnon peak is particularly evident along the $(\pi,\pi)$ direction, where the main magnon peak is sharp but notably asymmetric, which is not expected from the fitted model [64]. At lower $\mathbf{q}$ values and for the doped samples, the secondary magnon peak becomes less pronounced and may even vanish after background subtraction. The dispersive properties of this second peak can be observed for the $\mathbf{q}$ regions where its contribution has relevant weight to the fitting, close to the magnetic zone boundary, as shown in Fig. S11 of the Supplemental Material. This behavior supports the interpretation of the feature as a bimagnon, a bound state of two correlated spin flips on neighboring sites, with well-defined momentum and relatively long lifetime [65, 66, 67].

The peak energy of the magnetic excitation $\omega_{\text{peak}}$ takes into account the finite lifetime of the excitation by including the damping effect. We inspect it as a function of momentum, in Fig. 2(d), where it is possible to observe its energy dispersion as a function of transferred momentum for different dopings and directions. We found that the magnon dispersion is considerably softened for the Cr$8.5\%$ sample only, and equally softened in both directions, which contrasts with the Mn$8\%$ sample [68].

To further investigate this effect, we can also analyze the trend of the bare frequency, $\omega_{0}$, and the damping coefficient $\gamma$, as shown in Figs. 2(e) and (f). The $\omega_{0}$ $\mathbf{q}$-dispersion is nearly sample independent, and we compare it with the results of a linear wave theory, plotted with parameters obtained from inelastic neutron scattering experiments of the parent compound [62]. A good agreement is observed in the probed momentum region. The shadowed area represents the L dependency of the linear wave model, based on a Hamiltonian with a weak out-of-plane interaction component. The X/Y anisotropy of the detwinned sample is evidenced in this model, with our twinned sample dispersion lying in between the curves for the $\Gamma$X and $\Gamma$Y directions.

Along $(\pi,0)$, $\gamma$ is nearly $\mathbf{q}$-independent for the BFA and Cr$3.5\%$ samples, and becomes $\mathbf{q}$-dependent for Cr$8.5\%$, increasing monotonically with $\mathbf{q}$. It suggests shorter magnon lifetimes, likely due to competing Néel-type fluctuations derived from Cr. Along $(\pi,\pi)$, the damping varies non monotonically for BFA and, in particular, for Cr$3.5\%$. For Cr$8.5\%$, where the bimagnon peak can be neglected, the damping increases monotonically with $\mathbf{q}$, indicating that a genuine momentum-dependent $\gamma$ emerges for CrBFA with increasing Cr content. Our results align with studies linking Cr doping in CrBFA to weaker orthorhombicity and more localized (yet still itinerant) magnetism above $x>0.05$ [40] and theoretical modeling that increasing frustration may suppress the bimagnon peak in the FeSCs [67].

To understand the nature of the magnetically ordered state as a function of doping more comprehensively, we turn to ARPES measurements below and above $T_{\text{N}}$, as shown in Fig. 3. By inspecting the two different high-symmetry cut band maps, presented in the form of second derivatives of energy distribution curves (EDCs), we can access the band maximum points and see the electronic bands’ dispersion shape.

In the paramagnetic state, ARPES measurements of CrBFA reveal the increase in the hole pockets caused by Cr and the increasing electronic correlations, as characterized by the fractional scaling of the self-energy as a function of binding energy [24], a characteristic property of Hund’s metals [2, 30].

At $T=20$ K, the BFA sample shows the typical $d_{xz/yz}$ orbital splitting from stripe-type itinerant order [Fig. 3(a)], reflecting the itinerant spin character in the AFM phase. The same reconstruction is present for Mn substitutions up to $8.5\%$ along both $\Gamma$X and YZ cuts [17]. In contrast, for $8.5\%$ Cr substitution [Fig. 3(b)], no such reconstruction is observed. Instead, the bands crossing the Fermi surface around X (mainly $d_{xy}$) become noticeably flatter, especially in the ordered state.

The absence of the electronic structure reconstruction in the case of CrBFA implies a more localized nature for the magnetic ground state at this doping, as also observed in previous experiments [69, 40] and is consistent with the weaker orthorhombicity discussed above, since the $d_{xz}/d_{yz}$ splitting is intrinsically tied to the structural and magnetic phase transitions. Motivated by these observations, we developed and solved a theoretical XXZ model of localized spins to calculate the effect of the doping disorder on the magnetic excitations’ damping and dispersion, as detailed in the Supplemental Material S1.

Figure 4 shows theoretical results (dashed magenta line). Panels (a) and (b) present the dynamical structure factor, including intensities and linewidths (damping) as a function of Cr content $x$. We adopted $J_{1}=18.75$ meV to reproduce the magnon bandwidth at $x=0$. The dashed white rectangles mark experimentally accessible regions that are shown in more detail in the insets, wherein the color maps are rescaled to ease comparison with the experiment. Panel (c) shows the structure factor at $\mathbf{q}=\bar{Y}$ versus energy for different $x$, while panel (d) presents, for comparison, the experimental magnon, after elastic and phonon subtraction, for the same $\mathbf{q}$ point. The finite width of the theoretical curve at $x=0$ is due to the finite-size resolution of the calculations.

Our calculations explicitly include local fluctuations from Cr impurities and reproduce two key experimental observations: damping and softening of the excitations. The damping nearly doubles between $x=0\%$ and $8.5\%$, matching the experimental trend. The suppression of the energy scale of the excitations (the softening), however, is stronger in experiments than in the theory, likely due to the missing itinerant–localized dual character of the spin excitations. This indicates that charge doping also contributes to magnon softening. Panels (e) and (f) show the experimental $\mathbf{q}$ dispersion maps for BFA and Cr$8.5\%$ along the same high-symmetry direction as in the panels (a) and (b) insets. Within the experimental $\mathbf{q}$-range, theory and experiments agree with increasing spectral weight towards lower energies.

In addition to a good comparison with experiments, our calculations suggest that the magnon softening along the $(\pi,0)$ direction is related to a shift of Fe-derived spectral weight toward the $X/Y$ point with doping, as the magnetic excitations become more local, a $\mathbf{q}$-region inaccessible to RIXS. This theoretical insight may also help explain a key feature of the CrBFA and MnBFA phase diagrams: the absence of HTSC. Although Cr concentrates spectral weight near $(\pi,0)$/$(0,\pi)$, the associated excitations are highly incoherent (strongly damped), which is unfavorable for SC pairing [70, 71].

Other substitutional studies also indicate that local fluctuations from Cr moments can dominate over charge doping effects, as observed in Cr substitution of Ni-doped BFA [72, 73, 74] and P-substituted BFA [75, 76]. In CsFe2As2, disorder alone may even enhance correlations [77, 78].

Summary and conclusions: Our RIXS experiments investigated the evolution of the magnon dispersion along high-symmetry directions as a function of Cr content in CrBFA. Combined with theoretical modeling, our results suggest a scenario where, across the CrBFA $x$ vs. $T$ phase diagram, disorder dominates the damping and partially explains the softening of the excitations. It reveals that charge doping mainly contributes to suppressing the excitations. This places disorder, instead of charge doping, as the primary effect in CrBFA, similarly to the case of Mn substitution, where ARPES observes minimal hole doping.

This is unexpected given the increasing electronic correlation, changing hybridizations, and electronic filling as observed by ARPES experiments of CrBFA [24]. Our results show that disorder overrides these electronic effects in governing spin dynamics. We thus expect that our work may stimulate theoretical/experimental work on the effects of disorder on Hund’s metals.

S1 Theoretical model for magnetic excitations

In this section, we provide details about the theoretical model and numerical methods adopted in our investigation of the spectrum of magnetic excitations when disorder on the exchange interactions is gradually introduced. We describe the classical Hamiltonian and our approach to disorder, then we discuss how the dynamical spin structure factor (DSSF) at a given temperature can be numerically estimated by combining Monte-Carlo (MC) with semiclassical molecular dynamics (SMD).

S1.1 Magnetic model and disorder due to impurities

We consider the antiferromagnetic Heisenberg XXZ model [42, 43, 44, 45] on a square lattice for a minimal description of the magnetic excitations. We work on the classical limit, $|\vec{S}|\to\infty\hskip 1.13809pt$, and treat the localized magnetic moments (spins) in the lattice $\vec{S}(\mathbf{r}_{i})=S_{i}^{\mu}\hskip 1.13809pt\widehat{e}_{\mu}$ as unit vectors (with $N$ being the total number of sites, $i=1,2,\dots,N-1,N$ and $\mu=x,y,z$). The anisotropic Hamiltonian with exchange interactions up to next-nearest-neighbors is given by $H=\nicefrac{{1}}{{2}}\sum_{\hskip 1.13809pti,j,\mu}J^{\mu}_{ij}\,S^{\mu}_{i}S^{\mu}_{j}$, where $J^{\hskip 1.13809ptx,y}_{ij}=J_{ij}$ and $J^{\hskip 1.13809ptz}_{ij}=\lambda\hskip 1.13809ptJ_{ij}$ with $J_{ij}\geq 0$ and $0\leq\lambda<1$ so that the Hamiltonian favors the ordering of spins in the $xy$ plane (we consider $\lambda=0.95$). In the clean system of Fe-sites, $J_{ij}=J_{1}$ for nearest-neighbors bonds and $J_{ij}=J_{2}$ for next-nearest-neighbors bonds, otherwise $J_{ij}=0$. The system exhibits Néel order for $J_{2}/J_{1}<1/2$ and stripe order for $J_{2}/J_{1}>1/2$. In our study, we gradually introduce disorder in the exchange interactions for next-nearest-neighbors bonds $\langle\langle i,j\rangle\rangle$ by means of doping with Cr magnetic impurities as shown in Table S1. The values presented in the latter were motivated by experimental observations and estimates based on inelastic neutron scattering data [46]. As the RIXS resolution is limited in energy and momentum, this simple model is sufficient to capture the spectrum key features. For a given concentration $x$ of impurities, we generate disorder configurations by randomly selecting $\lfloor\hskip 1.13809ptx\,N\rfloor$ Fe-sites for substitution with Cr, each realization results in a list $D_{x}=\{\hskip 1.13809pti\hskip 1.13809pt\}$ mapping the impurity distribution in the lattice.

S1.2 Numerical calculation of the DSSF

The ensemble average of the energy-momentum resolved dynamical spin correlation function is the classical DSSF given by

[TABLE]

Here, $\mathbf{q}$ is a wave vector in reciprocal space, $\mu$ and $\nu$ are spin component indices, and $\mathbf{r}_{i}-\mathbf{r}_{j}=\mathbf{R}_{\hskip 1.13809ptij}$.

We follow a semiclassical approach, as discussed in previous works [47, 48, 49, 50, 51], to estimate this quantity using classical MC and SMD simulations: the former is employed to sample initial configurations of $N$ spins, $\{S^{\mu}_{i}(t=0)\}=\textrm{C}_{S}(0)$, and the latter to numerically integrate the equations of motion and obtain $\{S^{\mu}_{i}(t\geq 0)\}$. In our study, we consider a system of $N=L\times L$ spins, with $L=64$, at a fixed low-temperature value $T_{0}/J_{1}=10^{-2}$ and 4 different impurity concentrations $x\in\{0.00,0.05,0.10,0.15\}$. In our numerical approach for the $x>0$ cases, we generate a set of 100 disorder configurations $\{D_{x}\}$ and for each one a set $\{\textrm{C}_{S}(0)\}$ with 100 spin configurations is sampled from MC simulations. The MC implementation is based on the Heat Bath [52] and the Microcanonical [53] local update algorithms, where each MC iteration consists of 9 microcanonical steps followed by 1 heat bath step (here, a step represents a random sweep over the entire lattice). Independent simulations are performed for each $D_{x}$ configuration, after 50000 MC iterations for equilibration, we start drawing the $\textrm{C}_{S}(0)$ samples at every 100. To avoid twinned samples (i.e., with horizontal or vertical stripe modulation), we always initialize the MC code with a spin configuration exhibiting horizontal modulation. In Fig. S1, two $\textrm{C}_{S}(0)$ samples are depicted within a small region of the $64\times 64$ lattice. In (a), the concentration of impurities is small and the horizontal stripe ordering of the spins is barely altered. In (b), the concentration is high enough to induce a canting pattern that encompasses the whole lattice without disrupting the stripe order (moreover, small regions with higher local impurity concentration become common, enhancing the competition between Néel and stripe orders).

In a separate numerical implementation, we use the data sets collected to estimate the classical DSSF from Eq. (1) by means of the average

[TABLE]

Here, $\langle\hskip 1.13809pt\dots\rangle$ represents the ensemble average (i.e., the MC estimation based on the samples $\{\textrm{C}_{S}(0)\}\hskip 1.13809pt$) whereas $[\hskip 1.13809pt\dots]_{D_{x}}$ denotes the average over disorder configurations (absent for the clean case $x=0$). Each contribution to the combined average in Eq. (2) involves the Fourier transform in space and time $S^{\mu}_{\mathbf{q}}(\omega)$ of the temporal series $\{S^{\mu}_{i}(t_{n})\,|\,t_{n}=n\hskip 1.13809pt\delta t\hskip 1.13809pt,\,n=0,1,\dots,N_{t}\}$ for some spin configuration, where $N_{t}$ is the number of time steps within a finite time-interval with $t_{\textrm{max}}=N_{t}\hskip 1.13809pt\delta t$ representing the final evolution time reached under a fixed temporal discretization $\delta t$. The evolution equation can be derived from Heisenberg’s equation of motion for spin operators $d\hskip 1.13809ptS^{\mu}_{i}(t)/dt=\textrm{i}\left[H,S^{\mu}_{i}(t)\right]/\hbar$ in the classical limit [44], it describes the precession dynamics of each spin about the local field $h^{\mu}_{i}=\sum_{\hskip 1.13809ptj}J_{ij}^{\mu}\hskip 1.13809ptS_{j}^{\mu}$ produced by its neighbors:

[TABLE]

In SMD simulations, we employ the fourth-order Runge-Kutta (RK4) algorithm [54] combined with the energy correction (EC) algorithm from Ref. [55] to obtain $\{S^{\mu}_{i}(t_{n})\}$ given some $\textrm{C}_{S}(0)$ sample as initial condition. We use the Fast Fourier Transform (FFT) [56] algorithm to calculate $S^{\mu}_{\mathbf{q}}(\omega_{k})$, with $\omega_{k}=k\hskip 1.13809pt\delta\omega$ being discrete frequencies ($k=-N_{t}/2,-N_{t}/2+1,\dots N_{t}/2-1$). In order to obtain a suitable frequency resolution $\delta\omega=2\hskip 1.13809pt\omega_{\textrm{max}}/N_{t}$ over the entire frequency range while keeping $\delta t$ small for the RK4 + EC algorithm, we use $\omega_{\textrm{max}}=120$ and $N_{t}=4000$ which gives $\delta\omega=0.06$. The latter implies that $\delta t=\pi/\omega_{\textrm{max}}\approx 0.026$ and $t_{\textrm{max}}\approx 104.7$. With these parameters, our SMD code preserves very well the system energy $H(t)$ after each time step, with the final relative difference $\delta_{0}H=[H(t_{\textrm{max}})-H(0)]/H(0)$ of the order of $10^{-12}$ (without EC, the energy drifts away from $H(0)$ linearly and $\delta_{0}H$ is of order of $10^{-5}$).

We mention that the Fourier transform in time is applied to the input data $w(t_{n})\hskip 1.13809ptS^{\mu}_{\mathbf{q}}(t_{n})$, where the leading factor is a window function (usually a Gaussian envelope) that avoids numerical artifacts due to the FFT algorithm assumed periodicity. We use a custom window function that balances resolution and spectral leakage. It is designed to have a raised-cosine response with a very narrow transition from passband to stopband, and is approximated using a Kaiser FIR (finite-impulse-response) window/filter with appropriate parameters.

In order to compare our numerical results for the classical DSSF of Eq. (2) with experiments and theoretical approaches such as LSWT (linear spin wave theory), we follow the Refs. [57, 58, 51] and rescale $\bar{S}^{\hskip 1.13809pt\mu\nu}_{\textrm{\hskip 1.13809pt\tiny C}}(\mathbf{q},\omega)$ by the Bose-Einstein distribution $1/[1-\exp(-\beta E)]$ multiplied by $\beta E$, where $E=\hbar\omega$ is the excitation energy and $\beta=1/(k_{B}T)$ the inverse temperature. With this correction factor, we obtain the quantum-equivalent DSSF as (setting $\hbar=1$)

[TABLE]

The intensity plots (in arbitrary units) of the DSSF presented in the main text correspond to our results obtained from Eq. (4) for the longitudinal contribution $\mathcal{S}^{zz}(\mathbf{q},\omega)$. To enhance the visibility of spectral features across the energy-momentum regions of interest, we normalize the spectrum intensity and impose a high intensity cutoff at $20\%$ (when the region is the whole Brillouin zone, we instead impose a fixed low intensity cutoff prior to applying a logarithmic scale). For comparison with experiments, we use $J_{1}=18.75$ meV, which reproduces the $x=0$ magnon bandwidth, and it is compatible with inelastic neutron scattering results [46]. Additionally, we superimpose a dashed line indicating the LSWT energy dispersion that we derived for the Heisenberg XXZ model considered in this study. In Fig. S2, we show some results in two energy-momentum regions for two doping levels.

S2 Experimental Methods

Sample characterization: The samples were characterized by resistivity, specific heat, x-ray diffraction (XRD), and energy-dispersive x-ray spectroscopy (EDS) for their transition temperatures $T_{\text{N}}$, lattice parameters, and chemical content, respectively. The results were compared to composition vs. $T$ phase diagrams in literature [38, 39, 40] to benchmark the values of $x$.

RIXS experimental details: The samples were glued onto Cu sample holders using silver epoxy and cleaved in a vacuum of $1\times 10^{-8}$ Torr, using Al posts, at room temperature. The samples were then transferred to the measurement chamber, with a vacuum of $4\times 10^{-10}$ mbar. Experiments were performed in the magnetically ordered phase of BFA, with $T=15$ K. The scattering angle was fixed at $154^{\circ}$, for the maximum in-plane momentum transfer, and the sample angle was rotated so the momentum projection could be changed. The maximum momentum transfer was achieved for a grazing incidence geometry of $13^{\circ}$.

Linear horizontal (LH) and vertical (LV) polarized $X$-rays were used to probe the absorption spectrum (XAS) of each sample, and an energy-dependent RIXS map was done to study the incident energy dependence of the observed features. The energy of the absorption spectra maximum at $709.5$ eV and LH-polarized light with grazing incidence were chosen for the momentum-dependent studies. The RIXS experiments were carried out along the crystallographic directions $[110]$ and $[100]$. This corresponds, respectively, to iron sub-lattice first neighbors’ and second neighbors’ directions, also denoted as $(\pi,0)$ and $(\pi,\pi)$, or high-symmetry directions $\Gamma X$ and $\Gamma M$ of the reciprocal space, consistent with our previous ARPES results convention. For each RIXS spectrum, a carbon tape measurement was performed to determine precisely the energy resolution and zero energy transfer position.

ARPES experimental details: ARPES data were measured at the Bloch beamline of the Max IV synchrotron in Lund, Sweden. The total energy resolution was $\sim 10$ meV for incident photon energy of $76$ eV, and angular resolution of $0.1^{\circ}$. The LH polarized $\Gamma$X direction band map presented in this work is a high-statistic band measurement. The LV polarized YZ direction band map was reconstructed from Fermi Surface maps measurements. We present second derivative band maps for discussing band position, size and shape.

Fig. S3 shows the original intensity band map of the second derivative data presented in Fig. 3. For the BFA sample, it is evident that there is a band splitting as the sample is cooled, associated with the nematic orbital degeneracy lift.

S2.1 Low-energy RIXS background subtraction and fitting

All data presented in this study were normalized by the maximum fluorescence intensity before its subtraction, allowing a proper comparison of the quasi-elastic and magnon peak intensity as a function of momentum and Cr concentration.

The fluorescence background was fitted to a phenomenological model introduced by Yang et al. [59] as implemented already for the FeTe and BFA compounds [60, 61]. The model is illustrated in S4 and it separates the background into three regions with two crossover functions between them. The regions correspond to a linear behavior close to energy transfer $E\approx-1$ eV, shown as a yellow line, one exponential curve in the region $-4.5\lessapprox E\lessapprox-2$ eV, shown as a green line, and another exponential curve in the region $-6\lessapprox E\lessapprox-4.5$ eV, represented by a navy blue line. The complete model is depicted with the red line, which is a sum of the other three lines.

This model is described in equation 5, where the three curves are weighted by coefficients $\alpha$, $\beta$, and $\gamma$. In this model, the frequency $\omega$ corresponds to the transferred energy, which is defined as negative by $E=E_{d}-E_{i}$, where $E_{d}$ is the detected photon energy and $E_{i}=\mu h$ is the incident photon energy. The parameters $a$, $b$, and $c$ determine the exponential shape.

[TABLE]

The function that is responsible for the crossover between the curves is described in equation 6, where $\omega_{1,2}$ is the crossover position in transferred energy $E$ and $\Gamma_{1,2}$ determines the width of the crossover region.

[TABLE]

Using this function to subtract the fluorescence contribution allows us to have a longer zero line in energy loss for the resulting quasi-elastic plus magnon spectra, improving the fitting quality evaluation, such as residuals. To fit the quasi-elastic plus magnon resulting spectra, we use a model composed of two resolution-defined width Gaussians for the elastic line and phonon peak, one damped harmonic oscillator for the magnon peak, and a Gaussian for the secondary magnon peak.

In Figs. S5-S9 is possible to see all the resulting quasi-elastic spectra and their respective fittings, with all included contributions. The total fitting is the red line and the shaded regions are the fitting contributions, including the elastic line, phonon, main magnon, and secondary magnon peak, shown as pink, green, blue, and orange areas, respectively. We can observe the quality of the fitting when the background is properly subtracted, making the need for an additional magnetic peak evident. The result phonon is momentum and sample independent, with energy loss of 19(2) meV.

Figure S10 shows the trend of the fitted main magnon peak, normalized by the fitted intensity $I_{0}$. An offset is included to make the trend more evident, and the peak energy of the magnetic excitation $\omega_{\text{peak}}$ is marked by a red dot for each momentum. We can observe the dispersion and shape of the fitted magnon as a function of momentum, direction, and doping.

Figure S11 shows the dispersion of the bare frequency $\omega 0$ associated with the second magnon peak, with clear momentum dependence where its contribution is significant (orange peaks in Figs. S5-S9). The Cr$8.5\%$ data was omitted due to negligible relevance. At higher $\mathbf{q}$, the dispersion of this secondary peak follows approximately twice the main magnon energy (red line), supporting its identification as a bimagnon excitation.

Figure S12 shows the trends of the peak energy of the magnetic excitation $\omega_{\text{peak}}$, bare frequency $\omega_{0}$, and damping coefficient $\gamma$ for all the samples when the spectra are fitted using a single-magnon model. Under this constraint, the damping regime of Cr$3.5\%$ is closer to that of Cr$8.5\%$, which is qualitatively different from what is observed from the raw data. Another key difference is that the softening in the Cr$3.5\%$ sample becomes more anisotropic, resembling the behavior previously reported for Mn$8\%$.