Disorder- and Field-Induced Antiferromagnetism in Cuprate Superconductors

The underdoped high-Tc materials are characterized by a competition between Cooper pairing and antiferromagnetic (AF) order. Important differences between the superconducting (SC) state of these materials and conventional superconductors include the d-wave pairing symmetry and a remarkable magnetic response to nonmagnetic perturbations, whereby droplets of spin-density wave (SDW) order can form around impurities and the cores of vortices. In a simple picture, whenever SC is suppressed locally, SDW order is nucleated. Within a mean-field theory of d-wave SC in an applied magnetic field including disorder and Hubbard correlations, we show in fact that the creation of SDW order is not simply due to suppression of the SC order parameter, but rather due to a correlation-induced splitting of the electronic bound state created by the perturbation. Since the bound state exists because of the sign change of the order parameter along quasiparticle trajectories, the induced SDW order is a direct consequence of the d-wave symmetry. Furthermore the formation of anti-phase domain walls is important for obtaining the correct temperature dependence of the induced magnetism as measured by neutron diffraction.

A superconductor is characterized by a Bardeen-Cooper-Schrieffer (BCS) order parameter ∆ k (R), where R is the center-of-mass coordinate of a Cooper pair of electrons with momenta (k, −k). The bulk ground state of such a system is homogeneous, but a spatial perturbation which breaks pairs, e.g. a magnetic impurity, may cause the suppression of ∆ k (R) locally. What is revealed when SC is suppressed is the electronic phase in the absence of ∆ k , a normal Fermi liquid. Thus the low-energy excitations near magnetic impurities and in the vortex cores of conventional SC are essentially Landau quasiparticles trapped in bound states. The underdoped cuprates have been studied intensively in recent years in part because their proximity to the Mott insulator is thought to be responsible for many unusual properties, including possibly high-temperature SC itself. These systems are quite different from conventional SC, because when the pair amplitude is suppressed locally, e.g. by a vortex, a competing ordered state 1 stabilized by the proximity to the Mott state appears to emerge instead of a normal metal. This state is characterized at low temperatures T by static local SDW order with an ordering wave vector Q near (π, π), an order which is not present in the state above T c . This was reported first in elastic neutron scattering experiments 2 on La 2−x Sr x CuO 4 (LSCO), with a correlation length of several hundredÅ, but has been confirmed in other underdoped cuprate materials as well. [3][4][5][6] An enhancement of incommensurate static order was observed with increasing the applied magnetic field up to 14T. 2 Because the signal disappeared above T c , the magnetism was attributed to the vortices; indeed, scanning tunnelling microscopy (STM) measurements 7 on Bi 2 Sr 2 CaCu 2 O 8+δ (BSSCO) were able to directly image unusual charge order near the vortex cores which is almost certainly related to the field-induced SDW detected by neutron scattering.
Hints of magnetic ordering in the SC state had been detected earlier by µSR experiments in zero field 8,9 as a wedge-shaped extension of the "spin glass" phase into the SC dome of the temperature vs. doping phase diagram of cuprates (see Fig. 1a). Lake et al. 2 reported that an incommensurate magnetic order similar to the field-induced state was observed in zero field, too. But although it also vanished at T c , the ordered magnetic moment in zero field had a T dependence which was qualitatively different from the field-induced signal. The zero-field signal was attributed to disorder, but the relation between impurities and magnetic ordering remained unclear. Because strong magnetic fluctuations with similar wavevector are reported at low but nonzero energies in inelastic neutron scattering experiments on these materials, e.g. on optimally doped LSCO samples exhibiting no spin-glass phase in zero field, a b c d Figure 1: Phase diagrams and temperature dependence of magnetic order. a. Schematic temperature T vs. doping x phase diagram for cuprates. b. t/U dependence of the SC T c and the disorder-induced magnetic transition temperature T g for an impurity concentration n imp = 10% with potential strength V imp = 2.0t and a hole concentration x = 0.1. t is the nearest-neighbor hopping amplitude and U denotes the Hubbard on-site repulsion. c, d. Magnetic Fourier component |M (q)| 2 at the ordering wavevector integrated around (π, π) vs. temperature for (c) different interaction strengths U (at fixed V imp = 2.0t) and (d) different impurity potential strengths V imp (at fixed U = 2.5t).
it is frequently argued that impurities or vortices simply "pin" or "freeze" this fluctuating order. 10 Describing such a phenomenon theoretically at the microscopic level is difficult due to the inhomogeneity of the interacting system, but it is important if one wishes to explore situations with strong disorder, where the correlations may no longer reflect the intrinsic spin dynamics of the pure system. Such an approach was proposed by the current authors in a model calculation for an inhomogeneous d-wave SC with Hubbard-type correlations treated in mean field. 11 In this model a single impurity creates, at sufficiently large Hubbard interaction U and impurity potential strength V imp , a droplet of staggered magnetization with a size corresponding to the AF correlation length of the hypothetical pure system. [12][13][14] Such local impurity-induced magnetism has been studied extensively both theoretically and experimentally, and was recently reviewed in Ref. 15. As was shown in Ref. 11, when these droplets come close enough to interact, there is a tendency to form incommensurate, phase-coherent Néel domains whose size is sufficient to explain the observations by Lake et al. 2 in zero field. In addition, such a model explains the empirical observations that both increasing disorder 16 and underdoping 8 enhance the SDW order.
In this paper we investigate the origin of the "order by disorder" phenomenon described in Ref. 11, as well as the T evolution of the disordered magnetic state in applied magnetic field.
In the case of the field-induced SDW, an apparently very natural approach to the problem was developed by Demler et al., 17 who constructed a Ginzburg-Landau (GL) theory for competing SDW and SC order in a magnetic field. This order-parameter phenomenology describes correctly the reduction of condensation energy in the vortex phase of the pure SC, and leads to a phase diagram qualitatively consistent with experiments. 1 The GL approach, however, ignores the energies of the quasiparticles moving in the inhomogeneous state which also can crucially affect the competition between SDW and SC order in these materials at low T , as we show here.
In a d-wave SC without AF correlations, a bound state of an isolated vortex is found at zero energy 18 due to the sign change of the order parameter on quasiparticle trajectories through the vortex core. On the other hand, solutions of the Bogoliubov-de Gennes (BdG) equations describing coexisting d-wave SC and SDW order [19][20][21][22] show that this resonance is split by the formation of the SDW; that is, the system can lower the energy of the nearly bound quasiparticles by moving them below the Fermi energy. This finding is consistent with the STM experiments in the Abrikosov state of YBa 2 Cu 3 O 7−δ (YBCO) by Maggio-Aprile et al. 23 and of BSCCO by Pan et al., 24 which observed split peaks in the vortex cores. A similar bound state is associated with non-magnetic impurities such as Zn in BSCCO and has been also imaged by STM. 25 It is therefore important to explore the role of quasiparticle bound states and their coupling to the SDW order to identify the origin of both types of induced local AF in the SC state. A more detailed understanding of field-induced order is also highly relevant for the interpretation of the quantum oscillations observed in recent transport experiment in high magnetic fields. 26,27 These oscillations are possibly due to the formation of Fermi surface pockets as a consequence of SDW ordering and the concomitant reconstructed bandstructure. In addition, finding ways to understand the effects of disorder is crucial in order to reveal the intrinsic AF correlations present in the underdoped part of the cuprate phase diagram, not least because the magnetic fluctuations at higher energies may be responsible for SC itself.
The inhomogeneous mean-field theory presented here for electrons hopping on a square lattice with a d-wave pairing potential and subject to a Hubbard on-site repulsion U, reproduces the essential aspects of the field-induced spin-glass phase shown schematically in Fig.   1a. The primary purpose of this analysis is to model the inhomogeneous SC state and does not intend to describe the Mott transition to an insulating state at half-filling; therefore the doping dependence may not be directly compared to experiment. On the other hand, if doping is assumed to be correlated with the ratio of bandwidth to local Coulomb repulsion t/U in the model, a phase diagram very much like the one found in various cuprate materials is obtained, as shown in Figs. 1b-d. We consider this a reasonable qualitative approach, since changes in the Fermi surface of these materials reported by angle resolved photoemission spectroscopy (ARPES) over the "spin glass" doping range are small, 28 and it is therefore plausible that the primary effect on the electronic structure is due to the correlation induced band narrowing, as discussed in Ref. 11. The phase diagram we obtain in The problem studied here involves several length scales, in particular the inter-impurity separation, the inter-vortex separation, the SC coherence and the magnetic correlations lengths, which can be difficult to disentangle. We obtain results which reproduce well the qualitative aspects of the experiment by Lake et al., 2 but show that some features depend on nonuniversal aspects of disorder, in particular the process of domain wall nucleation.
While disorder-and magnetic-field induced SDW order both add to the ordered moment, the interference of disorder and magnetic-field effects is quantitatively significant. The domain wall formation proves responsible for the distinct T dependences of the field-and the disorder-induced magnetization. An intriguing aspect of the present theory is that it also includes a crossover from magnetic droplets to filamentary stripe-like structures in selected regimes of hole densities and impurity concentrations. The model therefore offers a concrete route to describe the physics of the pinning of stripe correlations in the SC state. This insight may prove relevant for many experiments in the underdoped cuprates which have been attributed to stripes.
The basis for our model analysis is the BCS pairing Hamiltonian for a d-wave SC with orbital coupling to an applied magnetic field, to which we add site-centered chemical disorder and a local Hubbard repulsion; the latter is treated in an unrestricted Hartree-Fock approximation: Here, c † iσ creates an electron on a square lattice site i with spin σ =↑, ↓. The hopping matrix elements between nearest and next-nearest neighbor sites i and j are denoted by t ij = t and t ij = t ′ , respectively. An electron moving in the magnetic field from site j to i acquires additionally the Peierls phase We start with a single non-magnetic impurity in a d-wave SC at U = 0. The fingerprint of the induced virtual bound state is a single near zero-energy peak in the local density of states (LDOS) (see Fig. 2a). In this situation there is no magnetization induced by the impurity. Increasing the on-site Coulomb repulsion beyond a critical value U c , a staggered magnetization emerges in the neighborhood of the impurity (Fig. 2b). The two-sublattice nature of the magnetic pattern, in conjunction with the spatial extent of the impurityinduced resonance leads to a spin-dependent splitting of the near-zero-energy peak in the LDOS; one resonant state with a selected spin direction is thereby shifted below and the other above the Fermi energy, as is most clearly seen at the nearest-neighbor sites of the impurity where the resonant state has the largest weight (see Fig. 2a). The spin-dependent The magnetization induced by an orbital magnetic field can be traced to the same microscopic origin as the impurity-induced magnetization. 20 Above a critical U c a staggered spin pattern is nucleated in the vortex cores with a spatial extent reaching beyond the size of a vortex core (see Fig. 2d), as observed in experiment. 2 For the parameter set chosen, the core radius estimated from the area where the order parameter is suppressed is about one lattice spacing. The LDOS in the vortex center reveals that the origin of the field-induced magnetization is tied to the spin-dependent splitting of the Andreev bound state in the vortex core (Fig. 2c). The conjecture that the field-induced magnetization indeed appears simultaneously with the peak splitting in the LDOS is explicitly verified in Fig. 3. For the unmagnetized vortex at T = 0.175t a single Andreev bound-state peak exists at zero energy.
With decreasing T the vortices nucleate a staggered spin pattern precisely at the T where the zero-energy peak in the LDOS splits. With further cooling the peak splitting grows, more spectral weight is shifted below the Fermi energy, and the magnetization is enhanced.
The natural next step is to consider a finite density of non-magnetic impurities in the presence of an external magnetic field and to compare it to zero-field results. Specifically for the modelling of LSCO, we assume in the following that the Sr ions are the primary source N i e iq·r σ z i (see the Supplementary information for details). The magnetic signal in the structure factor appears at the incommensurate wavevectors q = (π, π ± δ) and q = (π ± δ, π) (see Fig. 4a). The magnitude of the incommensurability δ however varies for distinct impurity configurations randomly selected for 22 × 22 lattice systems. For the weak magnetic signal in zero field the averaging over different impurity configurations is therefore imperative but computationally demanding. Applying an external magnetic field strongly enhances the magnetization and reinforces incommensurate peaks at unambiguously selected wavevectors which are robust against variations in the impurity configurations (See Fig. 4b).
Remarkably, the temperature dependence of the structure factor (see Fig. 4c) closely resembles the neutron scattering data on LSCO by Lake et al. 2 For the results shown in Fig.   4 we have chosen a parameter set where the staggered magnetization in zero field has its a b c Figure 4: Averaged magnetic structure factor for a many-impurity system. a,b. Intensity plot of the magnetic structure factor around (π, π) at T = 0.025t in zero magnetic field (a) and at finite field (b). The structure factor data were averaged over ten different impurity configurations. For the used system size of 22 × 22 lattice sites a magnetic flux of 2Φ 0 corresponds to a strong magnetic field with H = 59T . The impurity concentration n imp = x is fixed to 10% (V imp = 1.3t) and U = 2.9t > U c . c. T -dependence of the peak intensity integrated around (π, π) in zero-field and at finite field with the finite density of impurities n imp = x (blue and red curve, respecctively); for the data with Φ = 2Φ 0 and n imp = x the zero-field data were subtracted. For comparison also the structure factor in a clean system is included for the same magneticfield strength (purple curve). |M (q)| 2 (integrated) translates directly to the ordered spin moment squared in units of µ B per Cu 2+ . onset at a temperature T g indistinguishable from T c . This reflects a situation where upon cooling through T c the localized bound states inside the d-wave energy gap emerge and immediately split in the self-stabilizing staggered magnetic pattern. Towards lower T the magnetic structure factor rises in a markedly different way in zero and in finite field. While the field-induced part of the magnetic signal has a negative curvature, the zero-field magnetic structure factor increases approximately linearly upon cooling. The two mechanisms of a b c Figure 5: Anti-phase domain walls. a, b. Real-space image of the staggered magnetization at T = 0.025t induced by three non-magnetic impurities. In impurity configuration a anti-phase domain walls appear vertically. In configuration b the staggered magnetization induced by the three impurities adjusts to a uniform AF domain around them. c. Temperature dependence of the integrated magnetic structure factor for the impurity configurations a (blue curve) and b (red curve), respectively.
impurity-and field-induced SDW do not simply cooperate additively; the field-induced part of the magnetization is twice as large in the presense of impurities as compared to the fieldinduced SDW in the clean system (see Fig. 4c). The zero-field increase of the magnetization in the inhomogeneous SC state originates from the merging of AF patches nucleated by the individual impurities. Without impurities the increasing field-induced magnetization with cooling results from the growth of the magnetized regions around the well-separated vortex cores. Thus, both our zero-and finite-field results for the finite density of non-magnetic impurities n imp = x closely follow form of the T dependence of the the neutron scattering data for underdoped LSCO. 2 Still, due to computational restrictions we are not yet able to access the low magnetic-field strengths to allow for a direct comparison with experiment.
An important remaining question is why the T dependences of the magnetization are different in zero and in finite field. A hint is provided by the observation that the T dependence of the magnetic structure factor for the impurity-free field-induced magnetization and also for just two single impurities in zero field has a negative curvature. In both cases the induced staggered magnetization patterns around each impurity or each vortex, respectively, adjust their individual two-sublattice spin structures in phase and thereby avoid any domain walls. 29 For three nearby impurities, however, it proves already difficult to find a specific configuration where anti-phase domain walls are absent. In Figs. 5a,b we compare the staggered magnetization of two 3-impurity configurations with distinctly different domain-wall patterns. Remarkably, placing the three impurities on the same sublattice to form a rightangled equilateral triangle as in Fig. 5a generates a sequence of vertical anti-phase domain walls. If instead the impurities are configured in an acute equilateral triangle as in Fig. 5b a simply connected AF island forms around them. As Fig. 5c shows, with decreasing T the magnetic signal evolves differently for each impurity configuration. Intriguingly, |M(q)| 2 rises almost linearly for the configuration with vertical anti-phase domain walls while it has a negative curvature for the single domain island. These examples, and others we have investigated, suggest that the linear low-T rise of the magnetic signal for a finite density of impurities in zero field originates from the anti-phase domain walls which are always present in the randomly generated impurity configurations. For the field-induced magnetization it is the larger distance between the magnetized vortices which prevents the occurrence of domain walls and therefore alters the T dependence of |M(q)| 2 .
All the results presented above focused on static disorder-and field-induced SDW, but inelastic neutron scattering experiments have shown that in the SC state the spin excitations at finite energy have almost the same distribution of spectral weight in q as the frozen magnetic order. 30 For very low doping in the normal spin glass phase above T c , the neutron intensity pattern is rotated by 45 • and the connection to the spin correlations discussed here is less obvious. In fact, the utility of our model of choice is questionable for the description of the normal state where Fermi liquid concepts may not even be applicable. Nevertheless in the SC state we have provided a concrete foundation for the freezing of fluctuating spin correlations by disorder and magnetic field on the same footing; in particular the role of the quasiparticle bound states in the formation of the magnetic order has been highlighted.
The new picture that emerges is complementary to the global competition between SC a b Figure 6: Charge-density profiles. a. The same parameter set in zero field as in Fig. 4: U = 2.9t, V imp = 1.3t, doping x = 10% = n imp , and pairing interaction strength V d = 1.34t. b. U = 4.0t, V imp = 1.3t, doping x = 15%, n imp = 7.5%, and V d = 2.0t. In (a) and (b) the temperature is T = 0.025t. and AF phases in the sense that SC and disorder may significantly enhance SDW order in the underdoped regime. The d-wave pairing of the SC condensate is crucial for this generation of local magnetism, as we have shown. Support for this cooperative effect between SDW and SC comes not only from the onset of the elastic magnetic neutron signal at T c but also from Zn-substituted optimally doped LSCO. There it is found by µSR that 2% Zn induces a magnetic signal, but 3% Zn is found to eliminate it, but also destroys superconductivity 31 ; within the context of the current theory, this effect is understood not as a consequence of spin dilution 31 , but rather due to the destruction of the SC phase and thereby its ability to generate (or enhance) magnetic order.
Finally we show that a qualitatively different kind of inhomogeneous textures may also be stabilized within the present weak-coupling approach. Figure 6a shows the typical chargedensity profile in zero field for a parameter set used above to explore the onset of static AF.
As expected, at and near the impurity sites the electron density is reduced, and in these areas the local SDW patches nucleate. With increasing repulsion U and for larger hole densities exceeding the impurity cocncentration, the inhomogeneous spin and charge patterns change qualitatively, and the impurity-centered patches with reduced electron density evolve into hole-rich filamentary structures (see Fig. 6b). In this still SC state the filaments constitute snake-like paths through an SDW background with an average density of almost one electron per site. These textures provide a link to the study of disordered (quenched) stripes similar to those discussed recently within various GL models. 32 where the hopping amplitude between nearest neighbor and next-nearest neighbor sites i and j is described by t ij = t and t ij = t ′ = −0.4t, respectively. An orbtial magnetic field is represented by the Peierls phase factor ϕ ij = (π/Φ 0 ) r i r j A(r) · dr, while Φ 0 = hc/(2e) is the superconducting flux quantum and A(r) = B(0, x) is the vector potential of the magnetic field in the Landau gauge. The d-wave pairing potential is defined on two nearest neighbor sites i and j by where V d is the attractive pairing interaction strength, which we set to V d = 1.34t throughout the paper. We then define a gauge invariant d-wave order parameter on each lattice site i The chemical potential µ is adjusted to fix the average charge density n = 1 N i n i , while the electron number operator for spin σ at site i is given by n iσ = c † iσ c iσ , and the local charge-density operator by n i = n i↑ + n i↓ , respectively. S z i = 1 2 σ z i = 1 2 (n i↑ − n i↓ ) is the z-component of the spin-operator at site i. The Bogoliubov transformation diagonalizes the Hamiltonian in equation (A1), which thereby takes the form E 0 is the ground-state energy and γ † nσ creats an elementary fermionic Bogoliubov quasiparticle excitation with quantum number n, spin σ, and energy E nσ > 0. Calculation of the commutators of H from equation (A5) with the electron operators c iσ leads to a Schrödingerlike set of BdG equations with As we only search for solutions with positive E nσ , it is sufficient to solve the following single matrix equation This is because the solutions for E n > 0 are obviously identical to the solutions for the (positive) eigenvalues E n↑ of equation (A6) while for E n < 0 the following relation holds Since the solutions of equations (A6) and (A7) can be mapped on to those of the BdG equation (A9), we diagonalize equation (A9) to obtain the pairing potential ∆ ij , the charge density n i , and the local magnetization σ z i self-consistently from where f (E n ) = (1 + e En/k B T ) −1 is the Fermi distribution function and T is the temperature.
Sums over n run over positive and negative energies E n .
To maximize the size of the system for which equation (A9) can be diagonalized numerically, we take advantage of the magnetic translation symmetry of our model Hamiltonian (A1) by dividing the lattice into M x × M y identical supercells each with N x × N y sites (see Supp. Fig. 7). 20,35,36 We define the following magnetic translation operator 37 where R is the translation vector and T R translates any lattice vector r to the position r + R. Because [H, T R ] = 0, it is possible to block diagonalize the Hamiltonian H in equation (A1) using the eigenstates of T R . This reduces the eigenvalue problem (A9) of Applying the magnetic Bloch theorem block diagonalizes the BdG equations (A9), where k = 2π( nx MxNx , ny MyNy ), u nk (r i ) = u ink and v nk (r i ) = v ink . Thus we have to solve the following 2N x × 2N y matrix equation for each k where H ij and ∆ ij are only k dependent, if i and j belong to different supercells. Then the backmapping (see Supp. Fig. 8) leads to an additional phase for the u's and v's according to (A17), which is assigned to the matrix elements t ij (k) and ∆ ij (k). To make sure that two magnetic translations commute, we have to choose the magnetic field such that its flux through every supercell is a multiple of 2Φ 0 . 20,37 Hence 2Φ 0 provides a lower boundary for the magnetic flux threading each supercell, which corresponds for a supercell enclosing an area of e.g. 22a × 22a to a magnetic field of about 59 T (we assumed a typical value for the in-plane lattice constant a in the cuprates of about a = 3.8Å).
t' t R Figure 8: Hopping between supercells. A particle, which hops from the left supercell into the right supercell, is mapped back to the left supercell through the translation vector R. As a result the wave functions u and v obtain an additional phase given by the magnetic Bloch theorem (A17).
In order to make contact with neutron scattering experiments, we evaluate the magnetic structure factor S(q). In homogeneous systems it is defined as We approximate the spin-spin correlation function by the following factorization σ z i σ z 0 → σ z i σ z 0 . (A23) Because the system which we are interested in is in general inhomogeneous, we have to sum over all lattice sites. Hence we find the expression |M(q)| 2 = 1 N 2 ij σ z i σ z j e −iq·(r j −r i ) .
This approximation of the magnetic structure factor is identical to the Fourier transform of the magnetization squared.
In Supp.  Figure 9: Switching on an orbital magnetic field. a-c. Zero-field data. d-f. Finite-field data (Φ = 2Φ 0 ). a, d show the charge density n i , b, e the d-wave order parameter ∆ d i , and c, f the magnetization σ z i in real space. The same set of parameters is used here as in the rest of the paper, i.e. x = 10% = n imp , U = 2.9t, V imp = 1.3t. These data were obtained at the lowest temperature T = 0.025t we considered throughout the paper. Note the different scales in c and f.
Finally, for the parameters used here, the zero field case already contains impurity-induced antiferromagnetic order (see Supp. Fig. 9c), which is significantly enhanced by switching on a magnetic field. The magnetization peaks near the vortex cores, but due to the fact that strong type-II superconductors are penetrated by the field much beyond the cores, the magnetization is also enhanced in regions far away from the vortices, where the order parameter is nearly homogeneous. The SDW emerges due to the splitting of the Andreev bound state as explained in greater detail in the main body of the paper.