Nematic fluctuations in iron-oxychalcogenide Mott insulators

Nematic fluctuations occur in a wide range physical systems from biological molecules to cuprates and iron pnictide high-Tc superconductors. It is unclear whether nematicity in pnictides arises from electronic spin or orbital degrees of freedom. We studied the iron-based Mott insulators La2O2Fe2OM2M = (S, Se), which are structurally similar to pnictides. Nuclear magnetic resonance revealed a critical slowing down of nematic fluctuations and complementary Mössbauerr spectroscopy data showed a change of electrical field gradient. The neutron pair distribution function technique detected local C2 fluctuations while neutron diffraction indicates that global C4 symmetry is preserved. A geometrically frustrated Heisenberg model with biquadratic and single-ion anisotropic terms provides the interpretation of the low temperature magnetic fluctuations. The nematicity is not due to spontaneous orbital order, instead it is linked to geometrically frustrated magnetism based on orbital selectivity. This study highlights the interplay between orbital order and spin fluctuations in nematicity.


INTRODUCTION
Nematic phases occur in a variety of very different physical systems throughout nature. A nematic phase is formed when a discrete rotational symmetry is broken or reduced, while the translational symmetry is preserved. Nematicity has been studied mostly in the structural organization of soft matter such as liquid crystals [1][2][3] . In these materials it is easy to visualize the nematic state's defining characteristic: the presence of a preferred axis along which liquid crystal molecules are, on average, aligned. This spatial anisotropy is generally described by C 2 rotational symmetry. In the context of cuprate high-temperature superconductivity (HTSC) bond-nematicity was introduced as a competing order with HTSC itself. In recent years, nematicity has become the central focus of efforts to illuminate mechanism(s) of HTSC in iron-based superconductors (FeBS). By now, a variety of iron pnictides have been shown to exhibit electronic nematic states 4,5 and in numerous theoretical pictures, the nematic phase is thought to be intimately related to the superconducting mechanism in iron-based materials 6 The origin of the nematic phase in iron superconductors is heavily debated. Two scenarios are proposed 5 : One proposal suggests that anisotropic spin fluctuations are necessary precursors of the nematic ordering that has been experimentally observed 7 ; the other claims that ferro-orbital ordering involving d xz,yz orbitals is responsible for nematicity 8 . Complicating this debate is the fact that various Fe-based materials exhibit different magnetic order. Nematicity has been mainly studied in the 122 iron pnictides such as BaFe 2 As 2 , but recently nematic behavior was reported in the 1111 material LaAsFeO 1−x F x 9 . An important question emerges. Is nematicity a general phase behavior that exists in other Fe-based materials across the HTSC electronic phase diagram? For example, do Fe-based parent materials that are Mott insulators also exhibit nematicity? What role does electron correlation in multiorbital systems play in the formation of the interlinked nematic and spin density wave (SDW) phase and the competing superconducting phase 10 .
Increasing evidence points to an electronic mechanism of nematicity which would place the nematic order in the class of correlation-driven electronic instabilities, like superconductivity and density-wave transitions.
In this work, we report a combined experimental and theoretical investigation of magnetism intrinsic to the correlation-induced 11 Mott insulators La 2 O 2 Fe 2 O(S, Se) 2 (see Fig. 1). Neutron powder diffraction was employed to make a magnetostructural comparison of the materials. A short range nematic fluctuating behavior was revealed by pair distribution function methods. A critical slowing down of nematic fluctuations is observed in the nuclear magnetic resonance (NMR) spin relaxation rate; this result was complemented by Mössbauer spectroscopy data. In contrast to Fe-based superconductors, orbital nematicity does not play a role in the iron oxychalcogenides, since a d xz,yz orbital degeneracy is not present due to an alternating orientation of the FeM 4 O 2 octahedra. Therefore, a different mechanism must be responsible for the nematicity we observe in La 2 O 2 Fe 2 O(S, Se) 2 . Since these are Mott insulators, a strong-coupling view should be relevant.
Our magnetic neutron diffraction data reveal the establishment of antiferromagnetic (AFM) ordering and the similarity of magnetic behavior in these materials. The theoretical modeling explains the low-temperature nematic fluctuations as being related to quadrupolar fluctuations in a strongly frustrated magnet. This agrees well with NMR 1/T 1 data. The quadrupolar fluctuations are closely related to geometric frustration in the non-collinear spin structure, which itself arises as a consequence of strong orbital selectivity in the compounds. Thus, orbital-selective Mott correlations are important for understanding magnetism in the Mott-insulating states of the oxychalcogenides. Such selective-Mott correlations have recently been the focus of attention in the iron-selenides [12][13][14][15] , and the present study provides an example of the relevance of orbital-selective Mott physics in Fe-based materials.
The next section contains the results from different experimental techniques that were applied to La 2 O 2 Fe 2 O(S, Se) 2 . The experimental tools can be classified as either global or local structural and electronic probes. All work presented here was performed on powder samples. Neutron powder diffraction was performed in order to determine long-range magnetic ordering as well as structural ordering. The main text details the global magnetic ordering while Supplemental Material provides information about the global structural behavior. Nuclear pair distribution function (PDF), nuclear magentic resonance (NMR) and Mössbauer spectroscopy expose local structural details associated with shortlengths scale in proximity to specifically probed atoms. At the section's end, we present theoretical modeling that is helpful in interpreting our findings.

RESULTS AND DISCUSSION
Neutron diffraction Neutron powder diffraction experiments were performed on powder samples of La 2 O 2 Fe 2 O(S, Se) 2 using the C2 high-resolution diffractometer at the Canadian Nuclear Laboratories in Chalk River, Ontario. The C2 diffractometer is equipped with an 800 wire position sensitive detector covering a range of 80 degrees. Data were collected in the 2θ angular range from 5°to 117°using a Si (5 3 1) monochromator at a wavelength λ of 1.337 Å (see Supplementary Fig. 1 . We did not observe nuclear Bragg peak splittings, that would be indicative of a structural phase transition, in the La 2 O 2 Fe 2 O(S, Se) 2 diffraction data collected over the temperature range from 4 to 300 K. Thus, we did not detect a change of lattice symmetry over this temperature range.
We verified the consistency of the data with the non-collinear 2k magnetic structure previously obtained for La 2 O 2 Fe 2 OSe 2 16,17 . The high-temperature paramagnetic (PM) phase is compared to the 2-k antiferromagnetic phase in Fig. 2. The Sarah suite of programs 18 was used to analyze the representations and provide the magnetic basis vectors for refinement with Fullprof 19,20 . The magnetic cell is commensurate and is doubled, in both a and c, with respect to the structural cell. The ordering is associated with k 1 = (1/2, 0, 1/2) and k 2 = (0, 1/2, 1/2), and the single Fe site on (1/ 2, 0, 0) in the nuclear I4/mmm cell is described by two distinct orbits governing the two (1/2, 0, 0) and (0, 1/2, 0) Fe sites that are independent in the magnetically ordered state. However, we note that for powder samples the diffraction pattern is indistinguishable from the pattern assigned as the collinear AFM3 model in the literature 16,21 . Neutron diffraction shows that the magnetic structures of La 2 O 2 Fe 2 O(S, Se) 2 are similar with a distinction in the AFM onset temperatures T N . The magnetic Bragg peak intensity measured over a temperature range from 4 to 300 K, shown in Fig. 1, is proportional to the square of the magnetic order parameter. The data reveal the onset of AFM ordering at 90.1(9) ± 0.16 K and 107.2(6) ± 0.06 K for La 2 O 2 Fe 2 OSe 2 and La 2 O 2 Fe 2 OS 2 , respectively. The T N value for La 2 O 2 Fe 2 OSe 2 is in excellent agreement with values in the literature 11,16,21 .
Nuclear pair distribution function In order to gain insight on the short-scale geometrical atomic arrangment, PDF measurements on La 2 O 2 Fe 2 O(S, Se) 2 were performed using the nanoscale ordered materials diffractometer (NOMAD) 22 beamline at the spallation neutron source (SNS) at Oak Ridge National Laboratory (ORNL). PDF measurements are well-suited to provide details of the structural arrangement of atoms. Pair distribution function analysis involves obtaining the Fourier transform of the measured total scattering intensity S(Q) in order to obtain a real space representation of interatomic correlations. The pair distribution function is given as  Fig. 2 Schematic illustrations of magnetic phases. For spins (red arrows) located at Fe atom positions (blue circles), the arrangements are given for a the paramagnetic (PM) high-temperature phase. The low-temperature antiferromagnetic phase (b) is the non-collinear 2-k arrangement.
Q½SðQÞ À 1 sinðQrÞdQ, where Q is the scattering vector and r is the interatomic distance. PDF is a total scattering method i.e., meaning that both Bragg and diffuse scattering intensity data is simultaneously collected. Therefore, the technique is sensitive to deviations from the average structure 23 . Our PDF data were generated from the total scattering experiments using neutrons of a wavelength band from 0.1 to 2.9 Å. The maximum momentum transfer used for the Fourier transform was 31.4 Å −1 . The instrument parameters Q damp = 0.017 Å −1 and Q broad = 0.019 Å −1 were fixed for all of our fits and structural refinements of the data were performed using the PDFGUI 24 program and DIFFPY-CMI 25 suite.
We performed data refinements using an orthorhombic model of La 2 O 2 Fe 2 O(S, Se) 2 with a sliding fit range from (1.5-21.5) Å to (29.5-49.5) Å in 1 Å steps to investigate the temperature dependent structure on these various length scales (see Supplementary Fig. 2). This resulted in 29 fits for each temperature at which PDF data were collected. We have extracted the orthorhombicity parameter δ = (a − b)/(a + b) for both length scale and temperature series sequential refinements, where a, b are lattice parameters. Figure 3 shows the refined orthorhombicity δ extracted from fits of neutron PDF data. Structural models obtained from refining PDF data are strictly valid within the refined length scale 26 ; this specificity condition enables the probing of local structure on different length scales by varying the refinement range. By conducting such fittings, we observed higher orthorhombicity for short-range fitting (1.5-21.5) Å compared to long-range fitting (30.0-50.0) Å. This suggests that by tracking the orthorhombicity parameter we observe local scale distortions from tetragonal to orthorhombic structure. These spatially limited distortions represent fluctuating nematic order 27,28 over a broad range of temperatures. It should be emphasized that neutron powder diffraction, from the same samples, provides clear evidence that the average, long-range structures of La 2 O 2 Fe 2 O(S, Se) 2 remain tetragonal throughout the high and low temperature regimes 29 . However, the PDF results provide evidence of local symmetry breaking that results in nematic C 2 regions of fluctuating order 27,28,30,31 .

Mössbauer spectroscopy
In order to investigate the magnitude and the orientation of the ordered static Fe magnetic moment with respect to the electric field gradient (EFG) principal axis within the distorted FeO 2 S 4 octahedra, and to study the modification of the electronic surroundings via the electric quadrupole interaction, 57 Fe Mössbauer spectroscopy was performed on La 2 O 2 Fe 2 OS 2 and is discussed in comparison to similar work on La 2 O 2 Fe 2 OSe 2 17 . The 57 Fe nucleus probes the onsite magnetism of the Fe ion, therefore any significant changes in orbital and spin degrees of freedom should be reflected in the 57 Fe Mössbauer spectra.  The Mössbauer spectra were analyzed using the full static Hamiltonian with nuclear spin operators I z , I + = I x + iI y , and I − = I x − iI y . Here B is the hyperfine field at the 57 Fe site while Q, g I , and μ N indicate the nuclear quadrupole moment, g factor, and nuclear magneton, respectively. The polar angle θ and the azimuthal angle ϕ describe the orientation of the Fe hyperfine field B with respect to the EFG z-axis. The azimuthal angle ϕ was set to zero to avoid cross correlations. This is justified, because the asymmetry parameter η is small and the polar angle θ is close to zero for both systems. For both samples, the whole temperature series was fitted simultaneously to determine η, then η was fixed, and then the spectra were fitted individually. Within error bars, both samples In the magnetically ordered state a magnetic sextet is observed in La 2 O 2 Fe 2 O(S, Se) 2 data, indicating the presence of a magnetically long range ordered (LRO) state. In terms of the magnitude of the ordered moment (magnetic hyperfine field B (T)) and its temperature dependence, the nature of this LRO state is the same for both systems. The lines in Fig. 4b) indicate fits of the sub-lattice magnetization M(t) by using the equation MðtÞ / BðTÞ / ð1 À T=T N Þ β crit . The magnetic critical exponent β crit is deduced from these fits (see Table 1).
The 57 Fe Mössbauer spectroscopy provides clear evidence that in both systems the angle θ between the strongest EFG component (local z main axis), and the magnetic hyperfine field B hyp is equal to zero in the magnetically ordered regime, i.e., the magnetic hyperfine field B is oriented parallel to the local z-axis of the EFG principal axis system. The strongest EFG component is aligned along the O-Fe-O chains in agreement with local spin density approximation + U (LSDA + U) calculations 32 . This indicates that the Fe moments are oriented parallel to O-Fe-O chains. Therefore, both La 2 O 2 Fe 2 OS 2 and La 2 O 2 Fe 2 OSe 2 exhibit non-collinear magnetic order with a 90°angle between two different magnetic sublattices. This is consistent with three different spin structures as discussed for La 2 O 2 Fe 2 OSe 2 17 . The 139 La NMR experiments presented below decisively identify the 2-k structure shown in Fig. 2b as the actual magnetic structure in both systems. v (mm/s) 14 (7) 298 (4) The AFM transition temperature T N , low-temperature 57 Fe saturation hyperfine field B sat , and the critical exponent β crit are found to vary only in small ranges. Θ D is the Debye temperature.
B. Freelon et al.
The temperature dependence of the EFG is comparable in La 2 O 2 Fe 2 O(S, Se) 2 . Below ≈1.4 T N the principal component V zz decreases significantly in a near-linear fashion (Fig. 4c). Previously reported thermal expansion data of La 2 O 2 Fe 2 OSe 2 also showed an unusual change in this temperature range 21 . It is known 34 that thermal behavior can influence the values of the EFG tensors; however, another relevant possibility is that any degree of orbital electronic anisotropy can significantly modify the orbital occupancy. This would lead to a change in the valence contribution reflected in the EFG tensors. This particular scenario is well supported by the NMR spin-lattice relaxation rate data presented below.
Nuclear magnetic resonance Figure 5 shows the representative 139 La field sweep NMR spectra on La 2 O 2 Fe 2 OS 2 . Each spectrum consists of three pairs of satellites, and one split central transition. For lower symmetries, such as tetragonal systems, in the presence of finite EFG at the probed nucleus with nuclear spin I = 7/2, one would expect such NMR spectra. Upon lowering T, the spectral shape remains unchanged down to 109 K. At lower temperatures in the magnetically longrange ordered state the spectral shape changes considerably, because of the development of static internal magnetic hyperfine fields associated with the Fe ordering.
The field swept NMR spectra can be well described by diagonalizing the full nuclear Hamiltonian where B is the external magnetic field, z axis being the principle axis of EFG and γ 2π ¼ 6.0146 MHz T . B ⊥ (B ∥ ) are the hyperfine fields (internal fields) at the La site due to the Fe ordered moment. At temperatures T > T N NMR field sweep spectra can be adequately described without any additional internal fields (B ⊥ = 0 mT, B ∥ = 0 mT). In the magnetically ordered state two different local fields B 1 and B 2 at magnetically nonequivalent 139 La sites must be introduced in Eq. 2 to describe the observed spectra. The shape of the presented NMR spectra on La 2 O 2 Fe 2 OS 2 look very similar to the Se variant La 2 O 2 Fe 2 OSe 2 17 . We ascertain a non-collinear magnetic structure for La 2 O 2 Fe 2 OS 2 in agreement with the 57 Fe Mössbauer results. The presence of two different sub-spectra of equal intensity with local hyperfine fields parallel and perpendicular to the c-axis unambiguously confirm the 2-k magnetic structure to be realized in both systems 17 . The applied external field does not change the magnetic structure of the sublattices due to anisotropy fields much larger than 7 T. This was checked by 57 Fe Mössbauer spectroscopy, where a model based on pure superposition of internal and external field reproduces the measured spectra, similar to the NMR results.
There are no 139 La NMR spectral anomalies near or below 150 K in La 2 O 2 Fe 2 OS 2 nor near or below 125 K in La 2 O 2 Fe 2 OSe 2 17 . In contrast, the 57 Fe Mössbauer reveals an almost linear decrease of the principal component V zz of the EFG below 125 K and 150 K in La 2 O 2 Fe 2 OSe 2 and La 2 O 2 Fe 2 OS 2 , respectively. That no effect is observed in the 139 La NMR static parameters is not surprising because the local hyperfine interaction at the 139 La nucleus is not necessarily sensitive to the Fe electronic structure. Figure 6 shows the temperature dependence of La 2 O 2 -Fe 2 OS 2 139 (1/T 1 ) data collected at 5.0 T. The data are qualitatively similar to La 2 O 2 Fe 2 OSe 2 both in the paramagnetic and antiferromagnetic states 17   frequency (1/T 1 T~∑ q χ ″ (q, ω 0 )), where q is the wave vector and ω 0 is the Larmor frequency 35,36 . Therefore the NMR 139 (1/T 1 ) is very sensitive to fluctuations related to the orbital and spin degrees of freedom at the Fe site. The NMR 139 (1/T 1 ) data in Fig. 6 provide a clear indication of a continuous slowing down of the Fe spin or orbital degrees of freedom above the AFM ordered phase. Similar observations were also reported for LaFeAsO 9,37,38 where, just below the structural transition, 139 (T 1 ) starts to increase before it diverges at the antiferromagnetic ordering temperature. Both La 2 O 2 Fe 2 O(S, Se) 2 fit well into this framework albeit they have not been observed to undergo structural phase transitions. For T < T N , the spin-lattice relaxation rate 139 (1/T 1 ) is strongly reduced over a very short temperature interval. Below 87 K, for La 2 O 2 Fe 2 OS 2 139 (1/T 1 ) can be fitted by an activation gap-like behavior $ T 2 expðÀΔ=TÞ, with a gap Δ = 70 K, where as for the Se variant it is Δ = 55 K. This observation is consistent with the opening of a spin excitation gap as reported in a neutron scattering study 39 . For T < 30 K, we find a smooth crossover from gap-like behavior to a gapless regime, where 1/T 1 ≃ T 3 . Therefore, two distinct regimes are observed where 1/T 1 ≃ T 3 for T << Δ, while 1=T 1 $ T 2 expðÀΔ=TÞ for T > Δ.
The NMR 1/T 1 behavior for T < 30 K follow an T 3 power-law form which suggests the activation of an additional gapless spin fluctuation channel. We now address the origin of this behavior within the framework of an Heisenberg model with appreciable geometric frustration, an additional bi-quadratic term, and singleion anisotropy as it is appropriate for Fe-based systems.
Theory Motivated by earlier 40 local density approximation + dynamical mean-field theory (LDA + DMFT) calculations in which we estimated an S = 1 on each Fe site, we begin with the spin S = 1 bilinear-bi-quadratic model of the square lattice. This model serves mainly illustrative purposes (The current calculation provides an illustrative example as S = 2 for the high spin state for Fe atoms in M = S and Se as reported in ref. 16 ), and has been widely employed in context of magnetism and nematicity in Febased SCs 41,42 (While ref. 16 reported spin S = 2 to be consistent with inelastic neutron scattering data collected from La 2 O 2 -Fe 2 OSe 2 , we employed S = 1 in order to include the wellestablished bilinear bi-quadratic model plus the effect of magnetic frustration). The Hamiltonian is where S i is a spin-1 operator on site i. J ij = J 1a , J ij = J 1b such that J 1a is ferromagnetic (FM) and J 1b is AFM. J ij = J 2 for [ij] along the square-diagonal, and K is the coefficient of the biquadratic interaction for S = 1 along with D being the single-ion anisotropy term that fixes the orientation of the spin. It is known that this model admits various ordered phases, depending upon the relative values of the J ij and K for a given D 43 : a q = (π, 0) collinear antiferromagnetic (CAF) phase, a q = (π/2, π) AFM * phase, a q = (π, π) Néel AFM phase, and a q = (π, 0) anti-ferro-quadrupolar (AFQ) phase. Additionally, the strong geometrical frustration (supported by small value of the ordered moment relative to density functional theory estimates) inherent in the model, along with sizable biquadratic and single-ion terms can yield a non-collinear AF phase as observed in our systems. The FM and AFM couplings, along with sizable D, also favor D = 2 Ising-like spin fluctuations (as extracted from the critical exponent for M(T)). With this in mind, we have estimated the temperature dependence of the NMR relaxation rate, 1=T 1 ' lim ω!0 P q Im χðq;ωÞ ω by computing the local spin correlation function, χ ii ðωÞ ¼ R 1 0 dte iωt hT½S þ i ðtÞS À i ð0Þi using quantum-Monte Carlo (QMC) methods applied to the H [Eq.
Within the temperature range 30 < T < T N (~100 K), the experimental data for NMR 1/T 1 fits AT 2 exp(−Δ/T) to a very good approximation, with Δ estimated to be approximately 100 K. Moreover, at low temperatures T < 30 K, we also find that the 1/T 1 ≃ T 3 , indicating gapless spin excitations in the spin fluctuation spectrum. This is in line with earlier findings on the Se compound 17 . Our theoretical modeling is consistent with strong suppression of magnetic fluctuations below T N . The modeling also suggests, consistent with available inelastic neutron scattering data, a crossover to a gapped regime of spin fluctuations with a crossover scale O(50 K). Finally, at low T < 30 K, another spin relaxation channel opens up, indicating low-energy spin excitations below the gap.
These findings can be understood in by considering that a fieldtheoretic analysis for a collinear AF yields 1/T 1 ≃ T 3 for T << Δ, but 1/T 1 ≃ T 2 Δexp(−Δ/T) for T >> Δ. For an AFQ ordered state the same analysis yields T À1 1 ' T 3 e ÀΔ=T for T << Δ, but T À1 1 ' T 7 for T >> Δ. Both of these outcomes are inconsistent with our experimental findings and, at first glance, would rule out both collinear AFM and AFQ phases at low T in the oxychalcogenides. However, our results are not in conflict with the previously proposed non-collinear AFM phase in the iron oxychalcogenides 16,17 . Three-magnon processes can be ruled out, since no T 5dependence is observed in our experimental data (nor in the calculations) at any studied temperature 45 . In considering all of the presented data and our theoretical results, the clear scenario emerges in which the La 2 O 2 Fe 2 O(S, Se) 2 Mott insulators exhibit a non-collinear magnetic structure that might support the presence of AFQ phases. Furthermore, these materials each contain gapless excitations in the spin fluctuation spectrum. However, a detailed, full understanding of the low-temperature region T << Δ deviation Fig. 7 The nuclear magnetic resonance (NMR) spin relaxation rate (1/T 1 ). A large-lattice size Quantum Monte Carlo technique was used by employing the Algorithms and Libraries for Physics Simulations (ALPS) to estimate (1/T 1 ). By doing so, we computed the local transverse spin-spin correlation function for a S = 1 Heisenberg model with geometrical frustration, biquadratic exchange, and single-ion anisotropy as relevant for the Fe-based systems. The QMC results (green squares) were fitted with T 2 expðÀΔ=TÞ (magenta dashed line) and T 3 (blue line) in the high-T and low-T regions, respectively. The main panel indicates that a 1=T 1 ' T 2 expðÀΔ=TÞ form for T > 40 K, with an estimated spin gap of 98 K, smoothly evolves into 1/T 1 ≃ T 3 for T < 30 K. The low-T power law behavior emerges due to an additional two-magnon relaxation channel opening up at low T. Good accord between theory and experimental NMR data can be seen (see inset).
of the NMR 1/T 1 experimental data from a spin gapped regime needs to be developed. We mention that we are not able to rule out the possibility of impurities or defects causing 1/T 1 signal deviation from gap-like behavior at low temperatures. The presence of slight impurities was suggested to be the source of the low-temperature upturn in magnetic susceptibility data of La 2 O 2 Fe 2 O(S, Se) 2 29 . Were it possible to electron-dope the iron oxychalcogenides by chemical substitution or strain in order to modify their Mott insulating electronic structure, one would expect the small density of metallic carriers near the Mott transition to strongly scatter off these low-lying spin correlations. This would involve carriers propagating by emitting and absorbing two magnons in the usual view of a doped Mott antiferromagnet. The doped carriers would propagate by scattering off dynamical spin fluctuations that are encoded in the one-spin-flip fluctuation spectrum. We note that these processes might give rise to anomalous metallic features but they remain un-investigated. In addition, the orbital selectivity that arises from the multiorbital character of the iron oxychalcogenides will also generally lead to incoherent metallicity. How the interplay between these distinct electronic mechanisms might conspire to generate unconventional ordered state(s) in Mott systems remains a fascinating open issue. The current work should encourage efforts along these lines.
In summary, we have performed a combined theoretical and experimental study of the iron oxychalcogenides La 2 O 2 Fe 2 O(S, Se) 2 . A global structure analysis indicates the presence of a magnetic phase transition to a 2-k non-collinear AFM magnetic structure in both systems, while no average symmetry changes are detected. Striking signatures of local nematic fluctuations in the magnetic responses of these materials are observed as the AFM phase is approached from high temperatures. Nematic fluctuations are observed to occur by virtue of short lengthscale rotational symmetry breaking. The nematic fluctuation data is consistent with features observed with PDF and NMR exprimental techniques. Quadrupolar behavior is responsible for the critical slowing down of spin fluctuations was observed in the 139 (T 1 ) NMR data. These emergent quadrupolar fluctuations arise from sizable geometric frustration, which is itself a consequence of strong orbital selectivity of the compounds 46 . Our theoretical modeling shows that such orbital selectivity underpins numerous behaviors in the correlation-induced Mott-insulating states of the oxychalcogenides. Along with earlier DMFT work, the current results also establish the relevance of a strong correlation-based approach for iron oxychalcogenides. Finally, we note that nematicity has been invoked 47 in both experimental and theoretical contexts to describe order that competes with high-T c superconductivity in the cuprates and the iron pnictides. Therefore, the current work expands a current topic in iron-based high-temperature superconductivity while demonstrating the underlying relevance of structural nematicity to larger set of quantum materials.

Material preparation
Polycrystalline La 2 O 2 Fe 2 OM 2 M = (S, Se) samples were synthesized by a solid-state reaction method. The powders of La 2 O 3 (99.99%), Fe (99.99%), and Se or S (99.99%) were mixed according to the stoichiometric ratio, ground fully, then pressed into several slices. These slices were placed in an alumina crucible, sealed in an evacuated quartz tube, heated at 1000°C for 24 h, cooled to room temperature after the furnace was shut down. Then the sintered slices were fully ground, and the above process was repeated.

Neutron pair distribution function
The sample was contained in a standard 5.85 mm inner diameter vanadium can and mounted to vanadium tail liquid He cryostat. Data were taken between 2 and 300 K for a constant accelerator proton charge of 2.66 C, roughly equivalent to 40 min at full accelerator power. A diamond sample was measured for calibration and a vanadium sample was measured for normalization in the cryostat. Likewise, an empty container and empty cryostat measurement were used for background correction. 57 Mössbauer spectroscopy The 57 Fe Mössbauer spectroscopy was performed in transmission geometry on samples with an area density of 2.6 mg Fe/cm 2 prepared by evaporation of a methanol suspension of the powdered samples. The samples were mounted in a Cryovac He-4 bath cryostat. A commercial WissEl Mössbauerr spectrometer with a 57 Co-in-Rh source was operated in sinusoidal velocity mode. The transmitted 14.4 keV photons were detected using a KeTek SDD detector. NMR experiments 139 La NMR experiments in a temperature range of 5.0 ≤ T ≤ 170 K were performed by conventional pulsed NMR techniques with an optimized Ï/2 pulse condition. The NMR data were recorded in an 8 T superconducting magnet system with a 4 He variable temperature insert using a Tecmag LapNMR spectrometer. The polycrystalline powder sample was placed in a glass tube inside a Cu coil with a frequency of the resonant circuit 30 MHz. The spin-lattice relaxation rate was measured by using the saturation recovery method at the central transition in the fixed field. The data were well-fit to a single T 1 value. Both 57 Mössbauer spectroscopy and NMR experiments were performed on samples prepared from the same batch of crushed polycrystalline material.

Theoretical modeling
The quantum spin Hamiltonian was solved using the loop algorithm as implemented in ALPS. Convergence on the relevant observables were achieved by solving the quantum spin Hamiltonian on a 12 × 12 square lattice.

DATA AVAILABILITY
The data sets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.