Terahertz field control of in-plane orbital order in La0.5Sr1.5MnO4

In-plane anisotropic ground states are ubiquitous in correlated solids such as pnictides, cuprates and manganites. They can arise from doping Mott insulators and compete with phases such as superconductivity; however, their origins are debated. Strong coupling between lattice, charge, orbital and spin degrees of freedom results in simultaneous ordering of multiple parameters, masking the mechanism that drives the transition. Here we demonstrate that the orbital domains in a manganite can be oriented by the polarization of a pulsed THz light field. Through the application of a Hubbard model, we show that domain control can be achieved by enhancing the local Coulomb interactions, which drive domain reorientation. Our results highlight the key role played by the Coulomb interaction in the control and manipulation of orbital order in the manganites and demonstrate a new way to use THz to understand and manipulate anisotropic phases in a potentially broad range of correlated materials.

T he simultaneous ordering in multiple degrees of freedom is a general phenomenon in layered correlated solids. These ordered states are anisotropic and break the fourfold in-plane symmetry found in the high-temperature phase. In pnictides spin, charge, orbital and structural degrees of freedom lock together at the nematic transition 1 ; in the manganites charge, orbital, and structural degrees of freedom simultaneously order to give rise to orbital phases; similar transitions occur in cuprates 2 and nickelates 3 . Understanding which degree of freedom is primarily responsible for these transitions is a focus of condensed matter physics 1 and may help to explain why such phases are favoured over alternatives such as superconductivity. Furthermore, dynamic control of the anisotropy, which strongly influences the material properties 4 , could lead to novel devices based on correlated materials. La 0.5 Sr 1.5 MnO 4 (LSMO) is a prototypical manganite. The hightemperature state is a paramagnetic semiconductor. The crystal has a layered two dimensional structure of manganese planes as shown in Fig. 1a. Below the orbital ordering temperature, T OO ¼ 230 K, electrons in the twofold degenerate e g level of the Mn ions localize on alternating Mn sites. Simultaneously, the occupied e g orbitals align. This second-order phase transition breaks the in-plane symmetry, and two possible isoenergetic orbital domains can form. The domains consist of CE-type 5 zig-zag chains 6 oriented along one of two possible crystallographic directions, which give rise to anisotropic optical and electronic properties 7,8 . Below T N ¼ 110 K, antiferromagnetic order also emerges. The origin of orbital order remains debated, and electronic correlations [9][10][11] , structural distortions [12][13][14][15] , or spin interactions [16][17][18] have each been suggested as the primary factor.
In the following, we address this problem by focusing on the switching of one-domain type to the other under an applied THz field observed using optical birefringence. We show that the multisite Coulomb interaction provides the energy for domain reorientation, even when considering structural distortions. For the single-layer manganite LSMO studied here, field-assisted hopping of electrons from two highly occupied sites onto the same less-occupied site increases Coulomb repulsion and causes domain reorientation when the THz field is orthogonal to the CE chain direction. Our model also suggests that structural distortions create preferential electron hopping directions, which modify Coulomb interactions and may cause domain reorientation in other systems. As Coulomb interactions ultimately drive domain switching regardless of the field-coupling method, our results suggest that electronic corrections are a key factor in the formation of orbital order.

Results
Optically probing orbital order. We harness the anisotropic reflectance to measure orbital alignment. The polarization of an incident helium-neon laser (HeNe) beam is oriented to lie between the two zig-zag chain directions, A and B (Fig. 1b). Reflected light incident on a domain of Type A experiences a polarization rotation of opposite direction to that experienced by light incident on a domain of Type B. As the HeNe beam is much larger than the orbital domain size 19 , the spatially averaged reflected beam does not experience a net rotation when equal amounts of both domains are present (Fig. 1b). However, both domains have a non-zero projection onto the R > axis; thus the intensity of R > as a function of temperature allows us to observe the onset of orbital order, while R A -R B measures the preferential formation of one-domain type (Fig. 1c).
THz field-induced birefringence. Figure 2a shows the experimental set-up for detecting THz-induced domain alignment. Ten-ps-long THz pulses are generated at 13 MHz by a novel infrared free-electron laser (FEL) 20 . The polarization of the THz field is controlled by a mechanically rotating half-wave plate so that the polarization periodically lies along either domain direction A or B. Variations in the polarization rotation of a HeNe probe beam are detected at the rotation rate of the THz polarization using a lock-in amplifier. This enables the measurement of small changes in the probe polarization whilst maintaining a constant isotropic heat load on the sample (see Methods for details). We verified that the measured signal did not change when the rotation rate of the waveplate was changed from 1 to 4 Hz. Figure 2b shows the THz-induced anisotropy as a function of THz fluence for a sample temperature of 170 K, below T OO but above T N . The THz-induced anisotropy signal was observed to grow linearly with the field intensity for several different excitation wavelengths, and saturation occurred at the highest  intensities. By using the literature values for the anisotropy of a single domain 8 , we calculate that we could align B10% of the domains before saturation and heating effects dominate. Figure 2c shows that in this resonance-free region the induced anisotropy was weakly dependent on the THz wavelength. Due to the finite bandwidth of the FEL pulses, longer wavelengths have a longer pulse duration and lower peak fields. Thus, as the ability to align domains decreases slightly for longer wavelengths, we conclude that the peak field is more important for domain alignment than the total energy in the THz pulse, which remains constant.
To verify that the observed anisotropy signal was due to domain alignment by the field polarization, we performed several checks. Figure 2d shows the anisotropic signal as a function of the input HeNe probe polarization angle. As expected, when the probe polarization is rotated by 45°to align with an orbital domain axis the polarization state is not modulated by domain alignment, and a further rotation by 45°reverses the polarity of the signal. Figure 2e shows that a largely reduced signal was observed when an additional quarter-wave plate was inserted into the THz beam before the rotating half-wave plate to produce rotating circularly polarized light. The remaining signal is due to a small residual ellipticity in the polarization. This demonstrates the significance of the THz polarization direction and allows us to exclude pointing variations in the THz beam or thermal effects as the primary source of the anisotropy signal.
To understand these observations, we measured a detailed power and temperature dependence. Figure 3a shows the temperature and power dependence of power-normalized THzinduced net anisotropy. On cooling below T OO , the anisotropic signal increases rapidly for moderate THz powers. This occurs in step with the total orbital order shown in Fig. 1c. This rapid increase is due to the increasing correlation length of domains that occurs as the sample cools. On further cooling, the induced anisotropy saturates at B200 K and then starts to diminish, with the signal almost gone below T N . For higher THz powers the onset for anisotropy signal is shifted to lower temperatures. This shift can be identified as the result of DC heating, with 2 W of THz power inducing a heating of 85 K. The thermal origin of this shift is clearly seen in Fig. 3b, which shows the temperature dependence of induced anisotropy for different THz powers. At the lowest powers, identical power-normalized temperature dependencies are observed, demonstrating the linear and non-thermal nature of the signal. The high-power signal shows a large thermal shift. If this shift is compensated, the same power-normalized signal is observed. The shift in the rising edge can be used to determine the power dependence of the thermal contribution, which is found to scale as I 3 THz (white dotted line) and thus only contributes at the highest powers.
Measured THz field-coupling mechanism. The above results clearly demonstrate that the THz field can be used to align orbital domains, but the question remains as to how this is achieved. Thermodynamically, the energy of a dielectric material in an external field is proportional to e(o)E 2 where e(o) is the dielectric function at the frequency of the THz field. An anisotropy in the dielectric function at the THz driving frequency would create an energy difference between domains aligned parallel or perpendicular to the THz field. This energy difference is given as Although both cases have the potential for field-induced alignment, the switching mechanism will be very different. In the latter case the electric field does not couple to the order parameter directly and instead shifts the thermodynamic potential to favour one-domain orientation. Thermal fluctuations could move the material towards the more favourable orientation. Such a mechanism is unlikely to occur in our experiment as we use pulsed fields which are only 'on' for B10 ps. This makes this indirect process extremely unlikely as the majority of the time the sample is not in an applied field and the two domains are energetically degenerate.
Direct manipulation of the order parameter by THz fields has been achieved through a linear coupling of the field to resonant electric 21 or magnetic 22 dipoles. However, to date, domain control with THz fields has not been achieved experimentally. Furthermore, orbital ordering is not described by a macroscopic polarization or magnetization; thus, a linear coupling of the THz field to the order parameter cannot provide the driving force for domain reorientation in LSMO.

Hubbard model Hamiltonian.
To investigate how the THz field can still couple to orbital order, we consider a phenomenological 2D extended Hubbard model for the e g orbitals 1 j i ¼ d 3x 2 À r 2 and j i ¼ d 3y 2 À r 2 . The Hamiltonian, H, can be expressed as where U and V are the on-site and nearest neighbour Coulomb interaction terms respectively, n i,a denotes the occupation of the orbital a ¼ N, and n i ¼ n i,N þ n i, is the total occupancy of site i. The spatial anisotropy of the d-orbitals plays a key role in the manganites 6 , and for simplicity we take the limit that electrons can only hop between orbitals of the same character. The hopping between manganese ions decreases the energy of the electrons by t and occurs via a bridging oxygen 2p orbital. The hopping probability depends on y, the Mn-O-Mn bond angle, as t ¼ t 0 cos(p À y) 2 where t 0 is the value of the hopping for an undistorted, bond 23 .  Figure 4a shows the computed energy diagram for different orbital configurations as a function of the onsite Coulomb energy U. As we are studying the state of the system above the magnetic ordering temperature, the spin state is ignored. In this case, with to oU,V, the experimentally observed CE-type orbital ordering ground state is obtained from a fourth order perturbative expansion of the Hamiltonian, the details of which can be found in Supplementary Note 1 and are summarized in Supplementary  Table 1. As expected, the two types of CE-domains are energetically degenerate.
Previous studies on control of magnetic and orbital ordering in LSMO have focused on order melting from optical [24][25][26] or vibrationally resonant mid-IR pulses 27,28 . The former predominantly influences the charge system, whereas the latter modulates the lattice. Optical excitation triggers charge transfer between Mn sites, directly perturbing the orbital order. Resonant mid-infrared light couples to infrared active phonon displacement that directly modulates the Mn-O-Mn bond angle 29 , which modulates the electronic hopping energy, t. We now consider how these two mechanisms can break the degeneracy between the two lowest energy domains obtained from equation (1).
Competing coupling mechanisms. First, we consider the influence of the field directly on the charge degree of freedom. The THz photon energy is insufficient to trigger charge transfer excitations but instead polarizes the electronic orbitals along the field polarization direction. The effect of the THz field can then be treated quasi-statically in the length gauge as where r is the Mn-Mn bond length. We calculate the energy difference of the ground state when the field is applied along or perpendicular to the chain direction and find that the degeneracy of the orbital domains is split. The energy splitting was found to scale as Domains where the field is applied along the chain direction become more stable than those with the field perpendicular to the chains. The model recovers an energy splitting that scales as E 2 , matching the experimental observation. Details of the calculation are given in Supplementary Note 2. Second, we consider how the structure distortions induced by the ions in response to the electric field could induce domain alignment. As the oxygen and manganese ions are oppositely charged, they will move in opposing directions in response to the applied field, modifying the Mn-O-Mn bond angle and thus the hopping probability, t.
There is some controversy over the crystal structure of LSMO, with some groups reporting the I4/mmm structure above and below T OO (ref. 30) and others reporting of a structural transition to a Cmmm phase below T OO (ref. 31). In the high symmetry I4/ mmm phase the Mn-O-Mn bond angle is y ¼ 180°. As shown in Fig. 4b, the result of ionic motion in this phase uniformly decreases all bond angles in the same way. This renormalizes t to a lower value, but this does not split the degeneracy of the domain energies. Thus we do not expect to be able to control orbital domains in LSMO through structural distortions induced by the electric field in the I4/mmm phase or in any other material with linear bonds at equilibrium. However, in the Cmmm space group, the Mn-O-Mn bond angle is decreased to 176.72°at equilibrium. This distortion is, in general, larger in the cubic manganites. In this case, Fig. 4b shows that two bond angles are increased towards the ideal 180°angle and two bond angles are decreased under an applied field. Assuming small changes in bond angle, the change in hopping term becomes t ± ¼ t(1 ± aE tan y), where aE is the field-induced change in the bond angle y, which can be obtained from considering the response of a dipole in an electric field (see Supplementary Note 3). Equation (1) can be re-evaluated with the alternating hopping factors. The resulting domain splitting is given by The domain with lower energy is again the domain parallel to direction of the applied field, and the domain energy difference still scales with the intensity of the electric field.

Discussion
The domain-switching process that arises from these models are depicted in Fig. 5. In the field-driven case (a), the electric field pushes the charges between the Mn sites. When the field is aligned perpendicular to the domain chain direction, charges from two Mn 3 þ sites are forced on the same Mn 4 þ site. This increases the onsite Coulomb energy penalty and acts as the force to reorientate the domain. If the orbital flips, then the new domain structure is generated. In this case every Mn 3 þ charge is moved onto a separate Mn 4 þ site and thus does not experience the extra Coulomb repulsion. A very similar process occurs when the domains are aligned due to changes in the tolerance factor (Fig. 5b). Here when the field is aligned perpendicular to the domain direction, the hopping probability is improved such that charge from two Mn 3 þ sites is again moved onto the same Mn 4 þ . Like before, the Coulomb energy penalty acts as the force to reorient the domain by minimizing the energy penalty.
By considering the relative strength of the two processes we can determine the most likely driving mechanism in LSMO. Using values of U ¼ 5 eV, V ¼ 1 eV and t ¼ 0.25 eV, based on ref. 32, and the appropriate bond angles for LSMO in the Cmmm phase, we find that Dr t =Dr F j j% 1Â10 À 3 . Thus this simple model predicts that field-induced alignment is more likely for LSMO. Structural control is most effective when infrared displacements can induce large anisotropy in the hopping term. This occurs when the system starts in a state with large equilibrium bond angles, enabling larger angular changes for the same oxygen displacement. These displacements may be further enhanced by using fields resonant with the Mn-O-Mn bond. Thus the structural distortion mechanism may dominate in the cubic manganites, such as Pr 0.7 Ca 0.3 MnO 3 , where there are significantly larger initial Mn-O-Mn bond angles.
The model can also qualitatively explain the experimentally observed decrease in domain alignment efficiency when the system magnetically orders. Below T N the spins order ferromagnetically along the chains, but antiferromagnetically across the chains 6 . Although perpendicularly oriented domains are still energetically less favourable than parallel ones in the presence of the field, electrons that are transferred across the chains experience an additional energy penalty due to Hund's rule and domains can only switch if the localized t 2g spins also rotates. As the field does not couple to the spin degree of freedom directly, this process only happens due to random thermal fluctuations which are unlikely to occur during the short period of time in which the electric field is applied. Nearest-neighbour spin correlations appear well above T N and thus can start to impede domain alignment before long range spin order sets in (ref. 33).
Our results establish that pulsed THz fields can be used to control the anisotropy axis of manganites and represent the first demonstration of THz non-contact control over domains, moving beyond previous studies which examined melting of orbital order by THz fields. Our model shows the importance of Coulomb interactions as a driving mechanism for domain rotation in the manganites and strongly suggests that domain alignment is driven by field-induced hopping in manganites with low structural distortions. Combining the techniques presented here with X-ray based imaging experiments 34 would further verify the domain-switching process and enable observation of the domain-switching pathway and speed.
Similar experiments should also be applied to other anisotropic phases in the cuprates and pnictides to determine the relative strength of the Coulomb and structural mechanisms in these materials. Furthermore, control of orbital occupation is an emerging method for controlling electronic properties which has typically been achieved with strain or charge transfers from interface states 35 . Our results show that orbital control can be achieved dynamically with THz fields, opening new ways to control device properties and suggest orbital domains could be an alternative to magnetic data storage.

Methods
Birefringence measurement. Figure 1c shows the experimental set-up for detecting THz-induced domain alignment. A single crystal of LSMO was grown using the float-zone method, C-cut, optically polished, and mounted in a cryostat with quartz windows. Intense THz light and 632 nm light from a continuous-wave HeNe laser were incident on the LSMO sample at near-normal incidence (o0.5°). The polarization characteristics of the HeNe beam were analysed as a function of THz polarization.
Intense, tuneable, linearly polarized pulses of THz-radiation with Fourier limited pulse durations between 10 and 15 ps are generated at a 13 MHz repetition rate from the FELBE free-electron laser 20 . The THz light was focused to a diffraction limited spot onto the sample using an off axis parabolic mirror with through-hole to permit HeNe co-linearity. The THz polarization was rotated by a mechanical half-wave plate mounted in a motorized rotation stage set to rotate at 4 Hz. This causes a THz polarization rotation at 8 Hz. To check the effect of circularly polarized THz, a fixed quarter-wave plate was placed upstream of the rotating half-wave plate to give rotating circularly polarized THz. This modulation technique keeps the thermal load on the sample constant, while at the same time enabling the detection of the small domain alignment signal. Powers were measured with a calibrated power meter and the fluences reported were calculated by taking the pulse energy, obtained by dividing the power by the repetition rate of the FEL, and dividing by the spot size. THz light was attenuated upstream of the experiment using calibrated attenuators provided by the FEL facility. The THz spot sizes were measured by introducing a flip-mirror between the sample and the focusing off axis parabolic mirror and imaged on a pyroelectric camera (SPIRICON Pyrocam IIIHR).
Linearly polarized HeNe light was used to probe the effect of the THz on the orbital order. The polarization was chosen to lie between the two orbital directions such that there was a non-zero polarization projection into both orbital domains. The HeNe beam was focused to a spot much smaller than that of the THz on the sample, and the reflection was collected. The reflected beam passed through a beamsplitter to provide a reference of the total reflected intensity. The remaining HeNe light passed through a Wollaston prism to separate the light into Type A and Type B components and made incident on a balanced detector. The signal from the balanced detector was fed into a lock-in amplifier triggered at twice the frequency of the THz polarization rotation. This is appropriate as positive and negative THz electric fields produce the same force for domain rotation. The lock-in signal was measured over 10 s of acquisition time, and the mean value was used as the induced anisotropy. To normalize the wavelength dependence of the THz spot size, the induced anisotropy signal was normalized by the wavelength of the THz squared. The charge driven process. I When the field is applied perpendicular direction of the domain chain, the field causes charge from different Mn 3 þ sites to move towards the same unoccupied Mn 4 þ ion. This has an increased energy cost due to the Coulomb repulsion between the charges. II The energy penalty can cause electrons to change their orbital state so that each electron hops onto a different Mn 4 þ ion. III In the flipped state the orbital domain has been rotated. (b) The structural process. I When the field is applied along the chain direction, the Mn-O-Mn bond angle is increased and reduced asymmetrically. As a result hopping is increased such that electrons localized on different Mn 3 þ ions are more likely to hop onto the same Mn 4 þ site, again increasing the energy penalty in the field. II Again, the Coulomb energy penalty can be minimized when the orbital rotates and the new orbital pattern (III) is a rotated domain. Note that fields polarized along the crystallographic axis do not induce a domain reorientation for either mechanism.