Ubiquitous enhancement of nematic fluctuations across the phase diagram of iron based superconductors probed by the Nernst effect

The role of nematic fluctuations for unconventional superconductivity has been the subject of intense discussions for many years. In iron-based superconductors, the most established probe for electronic-nematic fluctuations, i.e. the elastoresistivity seems to imply that superconductivity is reinforced by electronic-nematic fluctuations, since the elastoresistivity amplitude peaks at or close to optimal Tc. However, on the over-doped side of the superconducting dome, the diminishing elastoresistivity suggests a negligible importance in the mechanism of superconductivity. Here we introduce the Nernst coefficient as a genuine probe for electronic nematic fluctuations, and we show that the amplitude of the Nernst coefficient tracks the superconducting dome of two prototype families of iron-based superconductors, namely Rh-doped BaFe2As2 and Co-doped LaFeAsO. Our data thus provide fresh evidence that in these systems, nematic fluctuations foster the superconductivity throughout the phase diagram.


INTRODUCTION
Unravelling the interaction of electronic orders in the phase diagram of copper-based 1,2 and iron-based materials 3 , such as stripe order 4 or nematicity 5-7 on the one hand and superconductivity on the other hand, became one of the most important tasks in understanding high-temperature superconductors. Most families of iron-based superconductors show a nematic phase in close vicinity to the superconducting transition, with nematic fluctuations being present in large areas of their electronic phase diagram. Therefore, nematic fluctuations have been widely proposed to be coupled to superconductivity and to be an essential ingredient to the pairing of electrons in these materials [8][9][10][11] . However, experimental evidence for nematic fluctuations in the overdoped superconducting regime is scarce and mostly limited to the underdoped regions. It remains thus elusive whether such fluctuations are an indispensable ingredient in the Cooper pairing or play the role of additionally enhancing the superconductivity, potentially driven by antiferromagnetic fluctuations.
The Nernst effect describes the occurrence of a transverse electric field E y perpendicular to a temperature gradient ∂ x T j j and perpendicular to an external magnetic field B z with the Nernst signal N given by N = E y /|∂ x T| (see Fig. 1a). In the linear regime with respect to B z one introduces the Nernst coefficient ν as ν = N/B z . The Nernst effect is therefore a transverse transport probe that combines thermal and charge excitations and is strongly influenced by electronic states near the Fermi level. Hence, it is expected that the Nernst effect is sensitive to fluctuations of electronic order parameters and of the Fermi surface [12][13][14][15][16] . Indeed, in pioneering experiments, the Nernst effect has successfully been employed for probing stripe/nematic phases and pertinent Fermi surface reconstructions in cuprate superconductors 1,17,18 . We will show below that the Nernst effect allows probing nematic fluctuations in iron-based superconductors and that it provides profound information on the interaction of nematicity and superconductivity.

RESULTS Theory
To obtain a conceptional insight on the influence of nematic fluctuations on the Nernst coefficient ν in iron-based superconductors, we consider a model for nematic fluctuations, where conduction electrons can occupy the d xz and d yz orbitals of the iron atoms on a square lattice. Such a minimal model captures basic features of the band structure of iron-based superconductors and allows for a ground state with nematic order 10,19 . The model Hamiltonian reads 10 , i;a n i;a À g 2 X i n i;xz À n i;yz À Á 2 ¼ H 0 þ H g ; (1) where a, b = xz, yz denote the orbital indices, σ = ↑, ↓ the spin index, and n i;a ¼ P σ c y iaσ c iaσ is the occupation of orbital a on lattice site i. The first term in Eq. (2) describes the kinetic energy of the conduction electrons. In order to start from a situation which is relevant to a variety of iron-based superconductors we take the hopping parameters t ab ij from ref. 9 , leading to the usual band structure of the iron-pnictides consisting of one hole pocket around the Γ point and electron pockets around X-and Y-points of the Brillouin zone. The chemical potential μ of the system is included in the second term. A variation of μ changes the relative size of the hole and electron pockets in momentum space, and, therefore, the electron filling. A variation of μ can be achieved experimentally by doping or applying pressure.
Most importantly, the third term in Eq. (2) accounts for the nematic fluctuations. Due to its quadratic form and negative sign, this term energetically favours (proportional to the coupling strength g) a difference in the local occupation between the d xz and d yz orbitals. Thus, if one orbital at lattice site i is preferentially occupied, occupation of the orbital in perpendicular direction is unfavoured. Together with the hopping term, such an effective electron-electron interaction captures the basic property of nematicity. We emphasise at this point that despite the local character of the nematic interaction, our theoretical treatment of the model goes beyond single-site fluctuations. Due to the presence of the itinerant hopping term, which combines not only different lattice sites but also different orbitals, our analysis naturally takes into account non-local higher order processes. Since the rotational symmetry is not broken by the Hamiltonian (2), the introduced model can be thought of being relevant for an iron-pnictide material that is close to a nematic instability.
Using this model we find (see Methods section for details) that the Nernst coefficient ν is sensitive to the occupation difference of the orbitals and is strongly enhanced by a finite nematic coupling g, which is shown in Fig. 1b. Taken into account the typical energy scales of phonons (10 meV) and the nearest-neighbour hopping matrix element (1 eV), g/t 1 can be expected to be in the order of 0.01. This enhancement of ν can be qualitatively explained by the process displayed in Fig. 1(c). Due to the nematic coupling electrons that are moving along a temperature gradient and are occupying two different orbitals at the same lattice site will demand an additional energy cost proportional to g, which is avoided by a further movement of one electron to a neighbouring lattice site. The presence of a perpendicular magnetic field B ! jjz leads to a coupling between xz and yz orbitals and therefore to an enhanced hopping perpendicular to the temperature gradient as well as the magnetic field.
Nernst effect in Ba(Fe 1−x Rh x ) 2 As 2 Inspired by the above considerations, we performed Nernst effect measurements (experimental setup schematically shown in Fig. 1a) in all important regions of the electron doped phase diagram on Ba(Fe 1−x Rh x ) 2 As 2 single crystals. Rh-doped BaFe 2 As 2 represents a canonical electron-doped iron-based superconductor with a phase diagram and superconducting properties which virtually are identical to those of Co-doped BaFe 2 As 2 20,21 . In particular, the spin density wave (SDW) ground state of the parent compound, which sets in at the magnetic/structural transitions at T N ≈ T S ≈ 135 K is gradually suppressed upon Rh-doping, in favour of a superconducting state at x ≳ 2.5%.
To take into account the effect of entropy which leads to an intrinsic linear temperature dependence, Fig. 2a displays the temperature dependence of the Nernst coefficient divided by temperature ν/T for different Rh-doping levels in Ba(Fe 1−x Rh x ) 2 As 2 on a semi-logarithmic scale. At most instances, ν/T is large and positive, and steadily increases by 1-2 orders of magnitude upon cooling. Further inspection reveals that this steady increase is always present whenever a sample is in the tetragonal paramagnetic, non-superconducting phase, i.e. at T > T S for x ≤ 4% and T > T c for x > 4% (see Fig. 2b). In contrast, in the magnetically Fig. 1 Measurement setup and theoretical insights. a Schematic Nernst setup for measuring the Nernst coefficient ν. A resistive chip heater is used to heat up the upper side of the sample (red). The bottom is coupled to a thermal bath of controlled temperature (blue), therefore a temperature gradient arises in the ab-plane. The magnetic field is applied in c-direction and the Nernst signal is measured perpendicular to the temperature gradient in the ab-plane. b Calculated Nernst coefficient ν as a function of the strength g of the nematic fluctuations. Starting from very small values ν increases dramatically. c Schematic picture of the mechanism leading to the effect in panel (b). An external temperature gradient ΔT generates a preferred hopping of conduction electrons (left panel). If the neighbouring lattice site is already occupied, an intermediate state with double occupancy is created (middle panel). Due to the Lorentz force of the perpendicularly oriented magnetic field B ! the enhanced energy from nematic fluctuations is preferably reduced by removing the electron in the perpendicular orbital (right panel) leading to enhanced transverse transport.
C. Wuttke et al.
ordered orthorhombic phase either large positive or negative contributions occur in the Nernst effect. In the case of nematic order, an interpretation of these low-temperature contributions is elusive since a strong in-plane anisotropy, uncontrolled twinning, and possible additional contributions from highly mobile Diraclike fermions need to be considered [22][23][24] . On the other hand, the strong additional increase, which is well recognisable near T c of the superconducting samples, is well-known to result from vortex motion and therefore is of no further interest in the present study.
While the Nernst effect is in principle sensitive to fluctuating Cooper pairs above the critical temperature and is capable of tracking the transition from vortex liquid to a phase-fluctuating superconducting regime 15,16 , it is rather the strong enhancement of the Nernst coefficient at high temperatures, far above any superconducting fluctuation regime, which attracts our attention. This enhancement is strongest in the parent compound of both families, where the ground state shows long-range antiferromagnetic order rather than superconductivity. In view of our theory and the known presence of nematic fluctuations, we expect the amplitude of these fluctuations to be directly encoded in the Nernst signal. Thus motivated by this, we have investigated the doping evolution of several isotherms of ν/T as shown in Fig. 2c. The data reveal a strong initial reduction of ν/T at low doping levels which tends to saturate upon reaching superconducting doping levels (x ≤ 4%). However, a significantly larger ν/T occurs at x = 6.1% before being strongly suppressed at x > 8%. This non-monotonic doping dependence is especially enhanced at high temperatures (T ≳ 200K). Remarkably, as seen in Fig. 2d (please note the colour coding of the amplitude for better visibility), ν/T tracks the superconducting dome and shows a maximum near the optimal doping level x ≈ 6%. These experimental findings, therefore, suggest a relationship between ν/T and the superconducting transition T c , whereas our theoretical finding points out the sensitivity of ν/T to nematic fluctuations. In a nutshell, this comparison suggests that electronic-nematic fluctuations play an important role for enhancing superconductivity.

Nernst effect in LaFe 1−x Co x AsO
To verify the universal character of our result across different families of iron-based superconductors, we performed Nernst measurements on high-quality single crystals of Co-doped LaFeAsO 25 . LaFeAsO is one of the most famous representatives of iron based systems, since superconductivity of this material class was first discovered in F-doped LaFeAsO 26 . Another way of inducing a superconducting ground state in LaFeAsO is replacing Fe with Co 27 . In contrast to BaFe 2 As 2 , LaFeAsO is showing separated magnetic and structural transitions. Both are gradually suppressed as a function of Co content x and vanish in a region x ≈ 4.5-5%, a superconducting dome develops above x ≈ 5% with an optimal doping level of x = 6% 25,[28][29][30] .
In analogy to the 122 system, two temperature regimes are excluded from the analysis: (i) due to the formation of twinned domains and a strong in-plane anisotropy the data in the low temperature nematic phase is elusive and (ii) ν/T in the superconducting phase as well as close to T c is governed by well-known contributions due to vortex motion and superconducting fluctuations 15 . Hence, we focus on the steady increase of ν/T by more than one order of magnitude upon lowering the temperature in the tetragonal phase. Figure 3a shows data of ν/T vs. T for all LaFe 1−x Co x AsO compositions on a semi-logarithmic scale. Similar to Ba(Fe 1−x Rh x ) 2 As 2 , ν/T displays a large positive value for all doping levels at room temperature, which increases monotonically upon cooling as long as a sample is in the tetragonal paramagnetic, non-superconducting phase. The temperature dependencies of ν/T in Ba(Fe 1−x Rh x ) 2 As 2 and LaFe 1−x Co x-AsO are showing strong similarities as well. Apart from a slight difference in magnitude, the features obtained by the temperature-dependent measurements seem to be an universal property of the Nernst effect in iron-based superconductors. In fact, similar levels of electron doping (Rh doping in BaFe 2 As 2 and Co doping LaFeAsO) produce remarkably similar ν/T vs. T curves (for example, in both systems the doping level x = 2.5% is the only composition where a negative Nernst signal is obtained). The extremely similar behaviour in Ba(Fe 1−x Rh x ) 2 As 2 and LaFe 1−x Co x-AsO, i.e. the evolution of ν/T as a function of temperature and  20 . T nem data is obtained from elastoresistivity measurements on the same samples (compare "Methods" section).
doping, seems to be of universal character and is therefore strongly suggesting that the same mechanism is driving the large Nernst coefficient due to the presence of nematic fluctuations.
Furthermore, the doping-dependent data of LaFe 1−x Co x AsO (Fig. 3c, d) reveal a non-monotonic doping dependence as well, albeit some differences. On one hand, the variation as a function of doping is smaller compared to Ba(Fe 1−x Rh x ) 2 As 2 , though an enhancement of ν/T above the superconducting dome is also visible in LaFe 1−x Co x AsO, which suggests that a correlation between the maximum of ν/T and T c is present as well. On the other hand, the data for LaFe 1−x Co x AsO even seem to indicate a slight double peak behaviour as a function of doping, which resembles recent reports of elastoresistivity, i.e., an alternative well-established probe for nematic fluctuations 28,31-33 , in LaFe 1−x-Co x AsO (measured on a superset of the crystals that were used in this work) 28 .
The existence of the same double peak structure in two distinct probes suggests that this is a robust feature of the nematic fluctuations in this system. Motivated by the analogy between Nernst and elastoresistivity in LaFe 1−x Co x AsO, we checked the consistency of the two techniques in Ba(Fe 1−x Rh x ) 2 As 2 . Interestingly we observed a mismatch, namely, the Nernst coefficient peaks at the optimal doping, as mentioned above, while elastoresistivity peaks in the underdoped region and does not track the superconducting dome (see Fig. 9 in the "Methods" section), in agreement with previous reports on Co-doped BaFe 2 As 2 31,32 .

Discussion
Clearly, this observation implies that elastoresistivity and Nernst effect probe similar, yet distinct aspects of the nematic fluctuations, where apparent differences are likely related to subtleties of the considered systems. This conclusion is corroborated by recent elasto-Seebeck and elasto-Nernst data which possess a somewhat different coupling to the nematic susceptibility than the elastoresistivity 24 . In this context, it is important to consider that all the elasto-transport properties, including elastoresistivity, crucially depend on the electron-lattice coupling, since, by definition, they represent the response of an electronic system to an external uniaxial strain. This is not the case for the pure Nernst effect, which is measured without stressing the lattice. In this sense, the Nernst coefficient can be sensitive to nematic fluctuations even far away from the actual structural transition, where the electron-lattice coupling is weaker. We further emphasise that, on a more general perspective, both scattering time and Fermi surface distortions and fluctuations thereof have to be taken into account to explain the transport nematic phenomenology. Since thermoelectric properties, especially the Nernst effect, are known to be extremely sensitive to the latter 12,13 it is natural to consider an orbital contribution, which eventually is connected to Fermi surface distortions, as a first candidate to theoretically approach the influence of nematicity on the Nernst coefficient in our basic model. We point out, as is discussed in ref. 24 , that an anisotropic scattering time can be expected to play a further important role in these systems.
We mention further, that, besides a finite nematic coupling, another important ingredient in determining the electronic properties of iron-based superconductors is their multi-band nature. From the orbital point of view, our minimal model, including a finite nematic coupling, arises within a two-band system, therefore it is intrinsically a multi-band effect. However, one could also argue that the evolution of Nernst coefficient as a function of doping could arise either from a doping-dependent charge-carrier compensation or a doping-dependent scattering rate. It is well-established that the ambipolar flow of electron-like and hole-like quasiparticles enhances the Nernst coefficient, unlike the Seebeck and the Hall effect, where the contributions of carriers with different signs tend to cancel out. This peculiar property of the Nernst coefficient makes its doping dependence predictable, namely a progressive departure from the multi-band compensation due to doping should lead to a monotonic decrease of the effect. In our case, this scenario seems to be insufficient, since we observed a pronounced non-monotonic doping-dependence of the Nernst coefficient. We could now suppose that the doping dependence of the Nernst coefficient  could be caused by a change of the scattering rate with the Rh-or the Co-content. However, within our model, we should introduce a non-monotonic change of the scattering time up to one order of magnitude to produce an effect comparable to the influence of a finite nematic coupling (see Fig. 7 in the "Methods" section). Even admitting different scattering rates for the electron-and hole-like bands, such a huge change of the scattering time would be unusual and incompatible with the monotonic doping dependence of the resistivity of electron-doped iron-based superconductors 34 . Hence, although one cannot exclude a priori that a doping-dependence of the charge carrier compensation and of the scattering time can play a role, a finite nematic coupling seems to be essential to reproduce our data.
Finally, we mention that antiferromagnetic fluctuations in principle can give rise to a Nernst signal, too 12,35,36 . However, it is well established by nuclear magnetic resonance (NMR) that in electron-doped BaFe 2 As 2 and LaFeAsO the low-energy spinfluctuations are monotonically suppressed towards optimal doping 30,37 . In contrast to this, our Nernst effect data clearly exhibit a non-monotonic doping dependence near optimal doping. Furthermore, the measured Nernst coefficient in the tetragonal phase always is linear in magnetic field, in contrast to reports where a non-linear in field Nernst response is interpreted as evidence of an enhancement due to spin fluctuations 35 . Thus, while we cannot exclude spin fluctuations to contribute to the enhanced Nernst effect at low doping levels, such can be ruled as the origin for the enhanced Nernst effect near optimal doping.
In conclusion, we established the Nernst effect as a sensitive probe for nematic fluctuations, and we show that the amplitude of the Nernst coefficient in the normal state tracks the superconducting dome in the electron-doped phase diagram of two representatives of major families of the iron-based superconductors. Our results, therefore, imply that nematic fluctuations are an indispensable ingredient for enhancing T c in the iron-based superconductors. Furthermore, our analysis suggests the Nernst effect as a principal technique to probe nematicity, and thus it should be considered for future experiments in order to unravel the mysteries of electronic order in iron-based superconductors.

METHODS Theory
Minimal model. We consider a two-band tight binding model, where conduction electrons can occupy the d xz and d yz orbitals of the iron atoms on a square lattice. Such a minimal model captures basic features of the band structure of iron-based superconductors and allows for a ground state with nematic order 10,19 . The model Hamiltonian reads 10 , Àμ P i;a n i;a À g where a, b = xz, yz denote the orbital indices, σ = ↑, ↓ the spin index, and n i;a ¼ P σ c y iaσ c iaσ is the occupation of orbital a on lattice site i. The first term in Eq. (2) describes the kinetic energy of the conduction electrons. In order to start from a situation which is relevant to a variety of iron-based superconductors we take the hopping parameters t ab ij from ref. 9 leading to the usual band structure of the iron-pnictides consisting of one hole pocket around the Γ point and electron pockets around X-and Y-points of the Brillouin zone. The chemical potential μ of the system is included in the second term. A variation of μ changes the relative size of the hole and electron pockets in momentum space and, therefore, also the electron filling. A variation of μ can be achieved experimentally for example, by doping or applying pressure.
Most importantly, the third term in Eq. (2) accounts for the nematic fluctuations. Due to its quadratic form and negative sign this term energetically favours (proportional to the coupling strength g) a difference in the local occupation between the d xz and d yz orbitals. Thus, if one orbital is particularly preferred at lattice site i the orbital referring to the perpendicular direction is at the same time avoided. Together with the hopping term such an effective electron-electron interaction captures the basic property of nematicity. Since the rotational symmetry is not broken by the Hamiltonian (2) the introduced model can be thought of being relevant for an ironpnictide material that is close to a nematic instability.
The model (2) has been studied theoretically by random phase approximation (RPA) 9 and determinant quantum Monte Carlo (DQMC) 10 . As is well known, the RPA works well in a weak-coupling regime. It predicts a variety of orders for Eq. (2), depending on the value of μ. In the range between μ/t 1 = 0.2 and μ/t 1 = 2.5, where t 1 is the nearest neighbour intra-orbital hopping matrix element, the RPA predicts the onset of uniform nematic order for g c /t 1 ≈ 1.7. For μ/ t 1 > 2.5, susceptibility peaks at wave vectors (0, π) and (π, 0) indicating stripe order have been found, while for μ/t 1 < 0.2 the susceptibility peaks appear at wave vector (π, π) and therefore predict antiferroorbital (AFO) (antiferroquadrupole) order.
Evaluation of the Nernst effect. The Nernst coefficient ν = E y /(−∇ x T) which is measured in this work can be generally expressed in terms of the thermoelectric tensor α ab and charge conductivity tensor σ ab as In the following, we consider a weak coupling regime and show that already a small value of g (in comparison to the parameters of the kinetic energy) is sufficient to strongly enhance the Nernst coefficient ν. As discussed above, for such a small coupling regime, the properties of the model Hamiltonian (2) are well-captured by the RPA, which is a form of a mean field theory. Motivated by this result, we proceed to study the influence of nematic fluctuations on the Nernst coefficient (3) as follows: At first, we apply the mean-field approximation to decouple the interaction term in (2) and to solve the model. Based on this decoupling, we calculate the transport coefficients in Eq. (3) separately. Thereby we consider a relatively small magnetic field B that produces the Nernst effect such that perturbation theory up to lowest order in B is sufficient to describe the experimental conditions. Using these approximations, we show that the significant properties of the measured Nernst effect are reproduced by the theoretical description, which then justifies the simplifications made to solve the model. The first step is to decouple the interaction term H g of the model Hamiltonian (2) using mean-field theory. After a Fourier transformation to the momentum space, we obtain at first the following expression for the kinetic energy part of the Hamiltonian, where ε k;x ¼ 2t 1 cos k x À 2t 2 cos k y þ 2t 3 cosðk x þ k y Þ þ cosðk x À k y Þ Â Ã (5) The hopping matrix elements t 1…4 are fixed according to ref. 9 to the following ratios with respect to t 1 , t 2 /t 1 = 1.5, t 3 /t 1 = 1.2, and t 4 /t 1 = 0.95. The dispersions ε k,x and ε k,y are shown in Fig. 4a for a specific momentum cut. From now on any energetic quantity will be given in units of t 1 . The mean-field decoupled interaction term of Hamiltonian (2) has the same operator structure as H 0 and reads in mean-field approximation, Àd xy ðc y kxσ c kyσ þ c y kyσ c kxσ Þ i À gNðn 2 x þ n 2 y À 2n 2 h Þ; (8) where N is the total number of k points in the Brillouin zone and n x ¼ hc y ixσ c ixσ i, n y ¼ hc y iyσ c iyσ i, and n h ¼ hc y Considering the transport tensor coefficients in Eq. (2) within a semiclassical Boltzmann dynamics in the relaxation time approximation, we obtain the following expressions 38 , k;a À μÞn ka ð1 À n ka Þ; (12) where τ is the scattering time, e the elementary charge, and n ka ¼ hc y kaσ c kaσ i the momentum-dependent occupation numbers of orbital 'a'. The numbers v xa = ∂ε k,a /∂k x and v ya = ∂ε k,a /∂k y denote the components of the Fermi velocity with respect to orbital 'a'. Eqs. (9)- (12) have to be evaluated explicitly in order to calculate the Nernst coefficient ν. Thereby, analytical expressions for the momentum derivatives of the dispersions ε k,a can be derived directly from Eqs. (5)- (6). Here, only the hopping parameters t 1…4 enter. Thus, the only property which is determined by nematic fluctuations (represented by the parameter g) is the momentum-dependence of the occupation numbers n ka .
Expressions for the occupation numbers n ka can be found easily by diagonalizing the mean-field Hamiltonian H MF . In the following, we quickly discuss the corresponding approach. Due to the bilinear operator structure of the parts H 0 and H g;MF the model Hamiltonian can be at first written in the form where by use of Eqs. (4) and (8) we have introduced A diagonalization of H MF can be achieved by using the new fermionic operators where the momentum-dependent coefficients u k and v k have to fulfil the following properties, Under these conditions, we find for the mean-field Hamiltonian a diagonal operator structure, and the Hamiltonian is now written in the 'band' basis where 'e' denotes the electron pocket and 'h' the hole pocket in the conventional notation of the iron-based superconductors. The quasiparticle dispersions read, where the dispersion functions E ðe;hÞ k now describe the actual band structure (for example, measurable by ARPES) consisting of two bands with the typical Fermi surface topology in iron-based superconductors. Explicit expressions for the momentum-dependence of E ðe;hÞ k can be derived from Eq. (22) by use of the Eqs. (14), (15) and Eqs. (5)-(7). In Fig. 4a, the band structure is plotted for the g = 0 case.
We now turn to the actual calculation of the expectation values hc y ka c ka i and hc y kx c ky i which enter all transport coefficients in Eqs. (9)- (12) and also the local expectation values n a ¼ ð1=NÞ P k hc y ka c ka i and n h ¼ ð1=NÞ P k hc y kx c ky i. Since the Hamiltonian H MF is diagonal in the formulation (21), the above expectation values can be easily calculated when the c-operators are replaced by the quasiparticle operators α and β using Eqs. (17) and (18). The remaining expectation values are bilinear in the quasiparticle operators and reduce to simple Fermi functions with respect to the dispersions in Eq. (22). We obtain where f(E) = 1/(1 + e βE ) is the Fermi function and u k , v k are given by the Eqs. (19) and (20).
Numerical results. We now show by numerical evaluation of the expectation values (23)-(25) that the Nernst effect becomes dramatically larger in the presence of nematic fluctuations. For this, we have calculated the transport coefficients (9)-(12) and the Nernst coefficient ν for given sets of the parameters μ, g, and T. We start with a discussion how ν generally behaves when the strength g of the nematic fluctuations is slowly switching from g = 0 to finite but still relatively small values. Figure 5a shows the calculated ν as a function of g for a fixed value of the chemical potential, μ = 1, and a fixed small Fig. 4 Modelled dispersion and occupation. Dispersion and occupation numbers of electrons in the d xz (blue solid lines) and in the d yz (green solid lines) orbital states plotted in momentum space along the k x direction. The chemical potential is fixed to μ = 1.2. a Difference between the orbital dispersions and the chemical potential, ε k,a − μ, calculated from Eqs. (5), (6). Since V t k vanishes along the k x -axis the dispersions ε k,a correspond in the absence of nematic fluctuations, i.e. g = 0, to two bands leading to the well-known hole-and electron-like Fermi surface sheets around (0, 0) and (π, 0). b Orbital occupation numbers (per spin direction) as a function of momentum along the k x -axis in the g = 0 case and a finite small temperature k B T = 0.02. Due to the absence of interactions, these expectation values correspond to Fermi functions. c Same quantities for g = 0.2. To highlight their momentum variations due to the nematic fluctuations, the occupation numbers are plotted between 0.9 and 1.03 for the d xz orbital and between 0 and 0.1 for the d yz orbital.
temperature value of k B T = 0.01. One can clearly see that ν increases strongly already in a relatively narrow g range, 0 < g < 0.2, before it reaches a saturation. We particularly note that already at g = 1/8 the absolute value of ν becomes three orders of magnitude larger than its initial value in the interaction-free case. In Fig. 5b, c, the transport coefficients which determine ν via Eq. (3) are plotted separately. It can be seen that both the Hall conductivity σ xy and also the anisotropic Peltier component α xy (blue lines) grow rapidly with g. Their increase is even stronger than that of the conventional longitudinal transport components σ xx and α xx (red lines). In the zero-g limit the Hall conductivity is much smaller than σ xx but due to the stronger increase it exceeds σ xx for g-values larger than g ≈ 0.2. This specific behaviour of the transport coefficients influences the Nernst coefficient ν significantly since it is determined by the ratio (3) where a summation of different tensor components enter in different ways in the nominator and denominator. It leads for ν to a change from strong increasing to rather saturating behaviour appearing around g ≈ 0.2 where σ xx and σ xy coincide.
The strong impact of nematic fluctuations to the anisotropic transport coefficients can be understood as follows. A closer inspection of the Eqs. (11), (12) shows that the lowest-order effect of nematic fluctuations is determined by the momentum-dependent occupation numbers n ka . All other quantities, including dispersions ε k,a and their derivatives are independent of g. In Fig. 4b, c, we show the momentum-dependence of the n ka for g = 0 (panel (b)) and g = 0.2 (panel (c)) where the Nernst coefficient has already significantly increased according to Fig. 5a. In the interaction-free case the occupation numbers are just the Fermi distribution functions separating fully occupied states (at T = 0), defined by the property ε k,a < μ, from the unoccupied states. The factor n ka (1 − n ka ) in the momentum summation in Eqs. (11), (12) is non-zero only in a narrow k-range around the corresponding Fermi momentum k ðaÞ F (compare panel (a)) where ε k,a − μ is of the order of the thermal broadening k B T. This limitation of the contributing terms in the momentum summation is wellknown and typical for the conventional thermal and electrical transport of conduction electrons in usual metals. It is, for example, the reason why the Nernst effect is usually very small.
Switching on the nematic fluctuations g changes the situation dramatically, as is seen in panel (c). Due to the coupling between the two different types of orbitals, there is now an increase of the d yz occupation (green solid line) around the Fermi momentum k ðxzÞ F of the other orbital, i. e. the d xz orbital. Note that the d yz dispersion ε k,yz is far above the Fermi level in this momentum region and thus would not be occupied in the decoupled case at low temperature. This non-zero occupation of d yz , which is solely due to the nematic coupling proportional to g, is relatively small but it is not bounded to some region around the Fermi level. Instead it affects the whole momentum range where the dispersion of the other orbital d xz is below the Fermi level. This leads to a macroscopic number of k-points where n ka (1 − n ka ) is non-zero, and thus the number of contributing terms in the momentum summations (11), (12) increases dramatically. Furthermore, as a consequence of the occupation n kyz , the number of d xz electrons, n kxz , (blue solid line) in the same k-state Hall conductivity σ xy (blue line) as a function of g. The Hall coefficient grows more rapidly than σ xx and eventually exceeds its absolute value at g/t 1 ≈ 0.2. This causes for ν a change from a strong increase followed by a saturation. c Peltier coefficients α xx (red line) and α xy (blue line) as a function of g. Also here, the anisotropic component is dominant in the presence of nematic fluctuations.  ≠ 0) is comparable with the observed behaviour of the Nernst effect as a function of doping as it has been found by our measurements. For comparison, we also show the result without taking into account nematic coupling (g = 0), which shows a very reduced amplitude. b Nernst coefficient as a function of temperature for different values of μ and fixed g = 0.1. As observed in the experiments, we find a μ-dependent low temperature maximum followed by a rapid decrease and saturation at high T. The maximum appears when the thermal energy k B T is of the same order as the experimentally found value k B T c for the superconducting T c . is reduced by n kyz . Finally, the occupation of both orbitals changes within the thermal broadening to 1 when k ðyzÞ F is reached. Note that the same discussion holds for the perpendicular direction in the k-space but with the opposite orbital characteristics.
In conclusion, the nematic coupling leads in the momentum range where for zero-coupling one orbital is occupied, and the other one is empty to a fluctuation-driven 'flow' from the occupied to the unoccupied orbital, which is oriented towards the transverse direction. Since this correlation affects also higher energy states, it has a dramatic effect on the anisotropic transport coefficients, which particularly describe the response to a flow in the perpendicular direction. Figure 1c shows the schematic picture of the responsible processes leading to this effect. An external field (for example, temperature gradient) applied in x-direction generates at first a hopping of conduction electrons in x direction. This can lead to an intermediate state with double occupancy, including both orbital types. Due to the nematic interaction in Eq. (2), this state has an order of g higher energy than the initial state. The enhanced energy can be reduced again by either further hopping in x-direction or (most importantly) removing the electron in the d yz orbital from the doubly occupied site. Due to its orbital character this electron preferably hops in y-direction where the hopping direction is pretended by the Lorentz force of the perpendicularly oriented magnetic field. The resulting preferred transport in y direction explains the high sensitivity of the Nernst effect towards nematic fluctuations.
In Fig. 6, we show the calculated Nernst coefficient ν as a function of the two further free parameters, the chemical potential μ and temperature T. One can clearly see from Fig. 6(a) that ν is mostly enhanced in a rather broad range around μ ≈ 1.8. Since μ is closely related to the electron filling, its variation can be considered as a simulation of doping the material with electrons or holes. Having this relationship in mind, we argue that the calculated maximum of ν is consistent with the measured behaviour of the Nernst effect for samples with different doping. The calculated influence of temperature is shown in Fig. 6b for three different μ-values. This result is also in qualitative agreement with our experiments. We find a μ-dependent maximum at rather low temperatures k B T ≈ 0.01 which corresponds to the temperature scale of the superconducting critical temperatures measured for this material. At higher temperatures, ν decreases again, followed by a saturation or even negative values in the high-temperature limit.
Finally, to study the effect of the scattering rate to our observed effect, we have artificially varied the scattering time τ in Eqs. (9)- (12) such that the calculated Nernst coefficient in Eq. (3) remains constant. The result is shown in Fig. 7. This corresponds to the change of ν, which would be necessary to overcompensate the effect of the enhancement of the Nernst coefficient by nematic fluctuations. It is seen that the change of τ must be extremely large, i.e. about one order of magnitude.

Crystal growth
Single crystals of LaFe 1−x Co x AsO were grown via Na-As liquid-assisted growth using diffusion-controlled abnormal grain growth from polycrystalline LaFeAsO by the solid state crystal growth method (SSCG) as described in ref. 25 . Single crystals of undoped BaFe 2 As 2 were grown by a high temperature solution growth using an excess of FeAs as solution (self-flux), (see, e.g., ref. 39 and references therein). Single crystals of Ba(Fe 1−x Rh x ) 2 As 2 were grown out of excess (Fe 1−x Rh x )As using the procedures outlined in refs. 40,41 , and 20 .
Crystals were analysed using several complementary techniques as scanning electron microscopy (SEM) with energy dispersive X-ray analysis (EDX), Laue backscattering, powder and/or single crystal X-ray diffraction, and SQUID magnetometry measurements.

Nernst effect measurements
The measurements were conducted in a home-built device in a helium bath cryostat under magnetic fields in the range ± 15 T. The single crystals were prepared in an inert atmosphere to prevent degradation. The thermal gradient in the crystals has been applied parallel to the ab-plane, using a resistive chip heater. The magnetic field along the c-axis and the Nernst voltage was measured orthogonal to the two former directions but again in the ab-plane, using a Keithley nanovoltmeter. The measured voltage has then been antisymmetrized with respect to magnetic field. All measured samples show a linear behaviour of the Nernst signal as a function of field (compare Fig. 8, which exemplary shows N(B) for two different doping levels of LaFe 1−x Co x AsO) in the tetragonal state, which is of particular interest for our analysis.

Elastoresistivity measurements
We measured the elastoresitivity of all the samples of the Ba(Fe 1−x Rh x ) 2 As 2 (with x = 0, 0.012, 0.025, 0.04, 0.061, 0.087, and 0.0107) series, in order to obtain the nematic phase diagram (Fig. 9b) discussed in the main text. For the measurement, we adopted the standard procedure described in refs. 32 and 28 . The samples have been glued to a commercial piezoelectric actuator, oriented in order to have the crystalline [110] direction aligned along the main straining-axis of the piezo stack. The applied strain has been measured with a strain gauge attached to the same piezo actuator. Figure 9a shows the temperature-dependence of the strain derivative of the resistivity dη/dϵ measured on a sample with x=0.025, as an example. Every point has been obtained at fixed temperature, by scanning the voltage applied to the piezo actuator from −30 to 100 V and detecting the corresponding variation in η and ϵ. Due to the tiny applied strain (ranging from 0.01 to 0.1%), we have operated in a linear regime of η vs ϵ. The nematic temperature T Nem has been obtained as a fitting parameter by interpolating the temperature-dependence of dη/dϵ in the fluctuating regime with a Curie-Weiss function dη/ dϵ = C 0 + C/(T − T Nem ) (red curve in Fig. 9a) 28,32 . C. Wuttke et al. Fig. 9 Elastoresistivity measurements. Elastoresistivity of the Ba(Fe 1−x Rh x ) 2 As 2 series. a Temperature-dependence of the strain (ϵ) derivative of the relative variation of the resistivity (η = Δρ/ρ) for a sample with x = 0.025. The red curves show the Curie-Weiss fit in the nematic fluctuation regime above the structural transition. b Phase diagram of the nematic susceptibility of the Ba(Fe 1−x Rh x ) 2 As 2 series. T C , T N , T S , and T Nem represent the superconducting, the magnetic, the structural, and the nematic critical temperature, respectively. The colour scale reproduces the magnitude of dη/dϵ, while the light yellow region corresponds to the superconducting dome.