Raman imaging of ferroelastically configurable Jahn–Teller domains in LaMnO3

The Jahn–Teller (JT) effect, through geometric deformation of molecules or local ionic lattices, lowers the overall energy of the system by removing electron degeneracy from partially occupied orbitals. Crystal symmetry lowered by JT distortion inevitably creates multiple variants of elastic and orbital-anisotropic states. Visualization and control of the domain/wall textures create a cornerstone to understand various correlated phenomena and explore wall properties. Here, we report the real-space observation of JT phonon and orbiton-related domains in a LaMnO3 thin film using confocal Raman spectromicroscopy. The characteristic symmetries of the JT-originated Raman modes allow us to detect and visualize the local population and orientation of the JT planes. Combined with a crystal structural analysis, we find that the formation of ferroelastic domains with W or W’ walls provides the basic framework for understanding JT domain textures. Furthermore, we demonstrate the JT domains can be manipulated by applying local external stress. Our findings provide a useful pathway for mechanically-tunable orbitronic applications.


INTRODUCTION
Transition metal oxides exhibit a plethora of functional properties and quantum phenomena, including high-temperature superconductivity, colossal magnetoresistance, and electrochemical catalytic activity via strong inter-couplings among spin, charge, and lattice degrees of freedom [1][2][3][4][5][6] . Double or superexchange interactions were proposed to account for the interplay of magnetism and transport properties 7,8 , and the basic principles strongly rely on the presence of charge ordering 9,10 , Jahn-Teller (JT) distortion [11][12][13] , and orbital ordering (OO) [14][15][16][17] . The JT distortion involves the partial removal of the d-electron degeneracy of a transition metal cation in a ligand field, resulting in a reduction in lattice symmetry. Since a system with orbital degeneracy is electronically unstable, the surrounding anions are deformed to remove the degeneracy. Cooperative manifestation of such local anion distortions leads to a long-range coherence, known as the cooperative JT distortion, which is necessarily accompanied by a long-range order of anisotropic orbitals, that is, OO.
The cooperative JT distortion is largely responsible for the complex relationship between the electronic and magnetic properties of the OO states [17][18][19] . For example, epitaxial deformation of La x Ca 1-x MnO 3 thin films promotes JT distortion, suppressing ferromagnetic properties and stabilizing the antiferromagnetic insulating phase 20,21 . A ferromagnetic phenomenon due to a self-doping mechanism was also observed at the polar-nonpolar interface of LaMnO 3 (LMO)/SrTiO 3 (STO) heterostructures 22,23 . Moreover, doped manganites undergo an inherent phase competition, causing the quantum ground states to fluctuate spatially. Accordingly, the appearance of domains, walls, and spatially non-uniform textures is an intriguing feature that strongly influences electrical and magnetic properties 9,24,25 . Similarly, the metal-insulator transition observed in LMO under high pressure is also understood based on a phase-separated configuration consisting of spatial regions of JT distorted and undistorted MnO 6 octahedra 26,27 . Theoretical analysis predicts a conduction enhancement at the anti-phase JT domain wall in an undoped manganite 28 .
To unravel the role of JT distortion and OO in strongly correlated functionalities and implement various proposed phenomena 11,25,28 , the observations of cooperative JT distortions and concomitant OO states are necessary. The fundamental properties have been successfully explored using resonant x-ray scattering [14][15][16] or through detection of characteristic phonon modes using Raman spectroscopy 26,29 . However, real-space access to the domain and wall textures of JT distortion along with nanoscale engineering and domain texture control is still difficult to achieve.
We employ confocal Raman spectroscopy, which is one of the far-reaching imaging systems capable of collecting high-quality multispectral images with a submicron lateral resolution and is widely used in a variety of research fields [30][31][32] . Its applicability is attributed to technological advances in high-efficiency laser sources, low-noise detectors, effective Rayleigh filters, and highthroughput optics (Fig. 1a) 33 . As a fast, convenient, and nondestructive technique, Raman spectroscopy has been used to directly observe the vibrational/rotational states of given crystal symmetries 26,29,34 and to confirm the elementary excitations of an order parameter such as magnon or orbiton [35][36][37] .
The perovskite LMO (space group Pnma) is a canonical system for OO harmonized with a strong JT distortion, by the hybridization of the Mn 3d and O 2p orbitals, wherein four electrons occupy Mn 3+ d-orbitals with a high spin configuration (t 2g 3 e g 1 ) 38 . Two-fold degenerate e g orbitals undergo the JT distortion (T JT~7 50 K) 39 , which lifts the degeneracy and lowers the electronic energy by single electron occupation (Fig. 1b). The OO state is established by cooperatively distorted oxygen ion cages connected to each other via corner sharing, resulting in the alternative occupation of the d 3x 2 Àr 2 /d 3y 2 Àr 2 orbitals in the orthorhombic ab plane 14,15 as schematically shown in Fig. 1b. A-type antiferromagnetic order is established via the super-exchange interaction between ordered inter-site e g electrons at T N~1 40 K 38 . In the bulk LMO with a D 2h point group symmetry, group theoretical analysis leads to 24 Raman active modes (7 A g + 5 B 1g + 7 B 2g + 5 B 3g ) among 60 Γ-point phonon modes 29 . According to lattice dynamical calculations 29 , it was shown that the phonon B 2g (A g ) mode near 610 cm -1 (~490 cm -1 ) corresponds to the symmetric (antisymmetric) stretching of basal oxygens of the MnO 6 as shown in Fig. 1c. The phonon mode with A g symmetry is associated with both the JT distortion and the octahedral rotation, but the phonon mode with B 2g symmetry originates purely from the cooperative JT distortion (Fig. 1c). The crucial effect of the JT distortion on the symmetric stretching mode can be identified from the symmetry analysis based on the character table of D 2h point group (see Fig. 1d and Supplementary Table. 1). The symmetry operation of C x 2 on a lattice with antiferrodistortive JT distortion results in a reversal of the stretching directions confirming the character (χ) of −1, but it becomes +1 in the case without the JT distortion.
In this paper, we implement real-space imaging of the cooperative-JT-distortion and OO domains in an epitaxial LMO thin film using polarized Raman spectroscopy. The interrelationships between dense striped JT domains, surface topographic features, and orthorhombic crystal structures are carefully investigated in an integrated manner. We discuss the JT domain Fig. 1 Polarized Raman spectra and maps of a 150-nm-thick LMO thin film. a Schematic of the experimental setup for the confocal Raman spectroscope with backscattering geometry. b Schematic diagrams of a JT distortion and the resultant energy shift of Mn e g levels, and the concomitant OO of singly occupied e g orbitals with the cooperative octahedral deformation within the ab plane. The ab JT planes are stacked along c axis and spaced from each other with a relatively smaller pseudocubic lattice parameter. c Raman spectrum measured in (yy) polarization geometry. In the inset, the anion vibrational patterns of two octahedra are depicted for A g and B 2g modes. d C x 2 symmetry operation (two-fold rotation around a-axis) on the vibrating LMO cells under the B 2g symmetry results in the opposite phase of the vibration direction, denoted by the character (χ) of −1. e AFM image. f Raman map of the A g mode at~490 cm −1 . g Raman map of the B 2g mode at 610 cm −1 . The areas marked in e-g with solid and dashed circles indicate the identical regions, respectively. The dark area in the lower right corner of the images is an intentionally burned area using a strong laser intensity. Scale bar indicates 2 μm. construction in terms of ferroelasticity and then demonstrate the elastic manipulation of the JT domains by applying local external stresses. This discovery provides a useful route for orbital engineering of complex oxides.

RESULTS AND DISCUSSION
Domain structure in a LaMnO 3 thin film We prepared a~150-nm-thick LMO thin film on a (001) C -oriented Nb-doped (0.5 wt%) STO substrate (where the subscript ' C ' stands for the crystal axes of the STO substrate in the cubic index) by pulsed laser deposition under a high-temperature growth condition that maintains OO even in thin films 15 . Figure 1c shows the polarized Raman spectrum that is measured in the (yy) polarization geometry (i.e., both incident and scattered lights have a linear polarization along y axis) and averaged over the whole region of a topographic image in Fig. 1e. Here, the x, y, and z directions are taken to be the directions of the crystallographic orientations of the substrate, i.e., [100] C , [010] C , and [001] C , respectively. The Raman spectrum shows peak features that have been observed in a previous report 29 , and the symmetries of the peaks can be assigned as aforementioned. Raman spectral peaks in the higher energy regime (>800 cm −1 ) have been suggested to be related to the elementary excitations of OO 35 , but the possible involvement of a lattice mode has not been completely ruled out [40][41][42][43][44] . Figure 1e-g display an atomic force microscope (AFM) image and Raman maps of A g (~490 cm −1 ) and B 2g (~610 cm −1 ) modes over an identical region. As clearly seen in Fig. 1e, the domain structures of the LMO film are composed of topographically flat regions and modulating regions along the {110} C directions.
Considering the lattice parameter relationship between the bulk LMO (a bulk = 3.990 Å, b bulk = 3.990 Å, c bulk = 3.846 Å in terms of a pseudocubic unit cell) 39 and the STO (a STO = b STO = c STO = 3.905 Å) substrate 45 , i.e., a bulk > a STO > c bulk , we presume that the ab JT planes (Fig. 1b) are most likely vertical to the substrate and a mosaic mixture of the vertical planes with four-variant orientations is responsible for the topographic feature through reducing misfit strain energy with the substrate.
The Raman map of the A g mode at~490 cm −1 (Fig. 1f) has a nearly uniform distribution. In contrast, the Raman map of the B 2g mode at~610 cm −1 (Fig. 1g) exhibits clear stripe domains with alternating higher and lower intensities. It is worth noting that the flat region (marked in the solid circle) has clear stripe domains in the Raman map of the B 2g mode, while the topographically modulated region (marked in the dashed circle) does not have distinctive stripes. We will analyze the Raman signals to clarify the correlation with the topographic feature in the following sections.

Angle-resolved Raman spectroscopy
Because of the possible multiple orientations of the JT planes (ab planes) in epitaxial thin films 15 , we first assess the symmetry of Raman intensity quantitatively concerning the possible geometries of the JT plane based on the group theory. Based on the symmetry analysis, the Raman tensor R of the B 2g and A g modes can be expressed in the light polarization space as and Here, each tensor element α ij (i, j = x, y, z) interrelates the incident light intensity of electric polarization (j) with the Raman scattered radiation of polarization (i). When we set the bases of the matrix as the principal axes of the orthorhombic LMO unit cell (note the subscript ' O ' represents an orthorhombic notation as described in Fig. 1b), the two Raman matrices can be purely diagonal or offdiagonal, and we approximate that the diagonal components of R O (A g ) have an identical value, i.e., α xx ≈ α yy ≈ α zz ≡ α considering the perovskite unit cell has a shape of pseudocube. Accordingly, the B 2g mode can play a decisive role in determining the possible orientations of the JT plane. Using Euler angles composed by elemental rotations to transform the lab frame (L) described in Fig.   1a into the orthorhombic frame (O), S L!O , we can find the Raman tensor in the lab frame for the given orientation of orbital ordering planes, and Raman intensity can be obtained using the following equation Here, e i and e j represent the polarization vector components of the scattered and incident light, respectively. Based on our previous report 15 , we consider the symmetry of the Raman modes for the three representative orientations of the JT plane as candidates. Figure 2 shows the JT planes parallel to the (100) C , (010) C , and (001) C planes ( Fig. 2a-c) and the resultant angle-resolved Raman intensities for two different polarization geometries ( Fig. 2d-g).

Fig. 3
Experimental results of the Raman intensities. a-f Polar plots of the Raman intensities of the B 2g mode normalized by that of the A g mode at three different positions, marked in the Raman maps above the corresponding graphs. The Raman maps were obtained at 610 cm −1 in the (xx) polarization configuration. The polar plots measured in the parallel (a-c) and cross (d-f) polarization geometries are displayed. The solid circles from 0°to 180°are experimental data, and the empty circles from 180°to 360°are the repetition of the data from 0°t o 180°in order to construct the polar diagram. g The fitting functions consist of two planes perpendicular to each other. The blue (red)-boxed terms indicate (010) C ((100) C )-oriented JT planes, as displayed in the schematics above the formula. h The definition of the director, which is parallel to the JT planes. Scale bar in the Raman map is 1 μm.
Y.-J. Kim et al. Here, the azimuthal angle (φ) in the polar plots spans from the substrate crystallographic [100] C direction to the polarization direction of the normal-incident light (Fig. 2). Figure 2d, e display the Raman intensities of the B 2g mode for the given orientations of the JT planes. If the JT planes are perpendicular to the substrate surface, the Raman intensities in the parallel polarization geometry become With two-fold symmetry, the maximum intensity is obtained when the light polarization is aligned along the JT planes. Meanwhile, the four-fold symmetry can be obtained when the JT planes are parallel to the substrate surface as In the cross polarization geometry (the polarization of incident light is perpendicular to that of scattered light), all of the cases have four-fold symmetry, but the intensity is significantly reduced when the JT plane is perpendicular to the substrate, i.e., Meanwhile, the Raman intensity I k ðA g Þ or I ? ðA g Þ of the A g mode shows no dependence of the azimuthal angle. I k ðA g Þ has a nonzero value in the parallel polarization geometry (Fig. 2f), however, I ? ðA g Þ vanishes in the cross-polarization geometry (Fig. 2g) regardless of the orientations of the JT planes. This is the reason why the Raman map in the B 2g mode has a clear color contrast, in contrast to the uniform A g mode map in Fig. 1f, g. Figure 3 and Supplementary Fig. 1 show the experimental results from the angle-resolved Raman spectroscopy of the B 2g mode for three representative positions, which are indicated by the Raman maps represented above the polar plots. For the quantitative comparison, the intensities of the B 2g mode were normalized by those of the A g mode. Figure 3a-c ( Fig. 3d-f) show the polar plots obtained in the parallel (cross) polarization geometry. The scattered data from 0°to 180°are experimental data (solid circles) and the data from 180°to 360°( empty circles) are a repetition of that from 0°to 180°, to provide a better comparison between the polar diagram and our fitting curves.
In the bright region of the Raman map, the intensity has a twofold symmetry, with a spindle shape where the intensity is maximized when the light polarization is aligned along the [100] C direction (Fig. 3a). In the dark region, the polar plot in Fig. 3b has a similar shape to the one in Fig. 3a but the graph is rotated by 90°. The two-fold symmetries suggest that the majority of the signals are from the (100) C or (010) C JT planes, respectively. Meanwhile, the area with the intermediate intensity has four-fold symmetry (Fig. 3c). The four-fold symmetry seems to come from the (001) C JT planes or from the composite structures of the (100) C and (010) C JT planes. The former case is excluded by the results of Fig. 3d-f, since the polar plots for all positions in the cross polarization case have the same symmetry and similar Raman intensities. Therefore, the JT planes of our LMO film prefer to be perpendicular to the substrate, and the symmetry of the polar plots can be attributed to the relative domain population of the (100) C and (010) C JT planes. Therefore, we set our fitting model for the polar plots based on the Eq. (5) with possible in-plane rotation φ of the JT vertical planes. f ¼ f jj ¼ a jj cos 4 ðφ À φÞ þ b jj sin 4 ðφ À φÞ for paralell polarization geometry f ? ¼ a ? sin 2 2ðφ À φÞ þ b ? sin 2 2ðφ À φÞ for cross polarization geometry: ( The first (second) term in Eq. (8) describes the azimuthalpolarization-angle dependence for the (010) C ((100) C ) oriented JT planes, as illustrated in Fig. 3g. Our model fits all of the experimental data well. Since the Raman intensity is proportional to the domain population 34 , the directional preference of the JT plane within a specified position can be determined. In this regard, we define a local director, which is a vector-like quantity with two equivalent opposite directions (headless vector), to understand the spatial distribution described in Fig. 3h. Such a director is widely used as an order parameter in nematic liquid crystals, where the local preferential order of the molecules affects the phase of the liquid crystals. In our system, the director describes the preferred direction of the JT planes, and the length of the director indicates the relative population of two perpendicular domains.

Real-space director field
To construct the real-space textures of the directors associated with the local orientation of the JT planes, we first established Raman maps over the same region for various polarization angles in the parallel polarization geometry. To align the multiple Raman images pixel by pixel, we used a modified image registration algorithm for optimizing the conversion matrices among images 46 . Based on our analysis method, introduced in the aforementioned section, the Raman maps of the B 2g mode were normalized by those of the A g mode, as shown in Supplementary  Fig. 2. We carried out fitting using Eq. (8) to find the coefficients per each position and mapped out the directors spatially.
We successfully visualized the distribution of the directors, thereby revealing detailed features on the striped domain structure, as shown in Fig. 4a. Figure 4b shows one representative normalized Raman map of the B 2g mode, measured in the (xx) polarization geometry, and Fig. 4c shows the corresponding AFM image. In addition, an R 2 map is also presented to assess the reliability of our fittings for constructing the director map in Fig.  4d, where there is little large deviation of R 2 from the ideal value +1. To convince ourselves, we also plot the fitting curves with the raw data of the polar plots at three representative points (Fig. 4e).
The widths of the stripe domains in Fig. 4a are clearly seen to bẽ 250 nm. The orientations of the directors are alternately aligned to the [100] C and [010] C axes. At the boundaries of the stripe domains, the magnitude of the director becomes vanishingly small, and we hardly detect an oblique orientation from the crystallographic axes. This is because the usual wall width of ferroelastic domains is less than tens of nanometers 47,48 . According to the Abbe diffraction limit, our spatial resolution is limited by d = 0.61 × λ/N.A = 0.61 × 532 nm/0.9~360 nm (here, λ and N.A. indicate the wavelength of the light and numerical aperture). The small magnitudes of the directors result from the spatial averaging of two equal-portion domains due to a narrow width of 150 nm, as estimated from the topographic images (see Fig. 4c and Supplementary Fig. 3).

High-energy regime
Until now we have discussed on the fundamental phonon modes associated with JT distortion and the corresponding domain structures. Meanwhile, the high-energy Raman modes (>800 cm −1 ) were proposed as collective excitations of the OO, i.e., orbitons, and their presence was measured experimentally 35 stimulating academic discussion on more detailed identities and possible origins, such as mixed modes with both orbital and phonon characters due to the phonon-orbiton interaction 40,41 , and multiphonon scattering activated by the Franck-Condon mechanism for a vibrationally excited self-trapped orbiton state 42,43 . A consensus of the origin has not been established, but the orbital characteristic is included in the proposed mechanisms. Thus, the observation of the orbiton-related modes and their correlation with the underlying phonon modes provide useful information for understanding the domain structure of the local OO beyond the JT distortion.
With our analysis of the phonon modes in mind, we investigate the real-space distribution of the orbiton-related Raman peak intensities, which are involved in a local OO in our system. Figure  5a shows the zoomed-in Raman spectra for the high energy regime, normalized by the intensity of the fundamental A g phonon mode. The inset figure is the Raman map for the B 2g mode measured in (yy) polarization geometry. Three representative positions were considered, as marked in the inset of Fig. 5a. Two orbiton-related peaks were clearly observed. The labeled A and B peaks show different tendencies in the intensity variation depending on the spatial position, suggesting the A and B peaks have A g and B 2g symmetry, respectively. It is supported by the fact that the real-space intensity modulation of the orbiton-related B peak coincides with that of the fundamental phononic B 2g mode, as compared in Fig. 5b. Since the magnetic structure of LMO at low temperatures is strongly coupled to the orbital configuration via the superexchange interaction 38 , the characterization of orbital-related domains offers a useful insight into the magnetic domain texture of LMO.

Crystallographic structure
To elucidate the microscopic origin of the stripe domains, we characterized the crystal structure using x-ray θ-2θ scans and reciprocal space maps (RSMs). Figure 6a displays the θ-2θ scan exhibiting the (00 L) PC diffraction peaks of the LMO film without noticeable impurity peaks (hereafter, ' PC ' stands for the pseudocubic indices of LMO film). From the θ-2θ scan, the outof-plane lattice parameter of the LMO film was estimated to be 3.98(6) Å, which is almost the same as the bulk in-plane lattice parameters in terms of a pseudocubic unit cell 39 . This implies that the orientation of the JT plane in our film is perpendicular to the substrate surface 15   In Fig. 6b, the structural modulation is reflected in a ω-rocking curve measured along [110] C at the LMO (001) PC , (002) PC , and (003) PC reflections. In addition to strong central peaks, two satellite peaks were observed at the identical off-axis ω-angles of ±1.1(7) (°) regardless of the out-of-plane reciprocal positions, which suggests the mosaic rotation of the unit cell rather than structural periodicity along the [110] C direction.
To clarify the crystal structure including in-plane lattice parameters, we obtained RSMs for two asymmetric reciprocal positions around the (103) PC and (223) PC peaks (Fig. 6c, d). To interpret the RSMs, we first considered a tilted pseudocubic unit cell (M C -type monoclinic unit cell of which the c-axis tilts toward <100> PC direction 49 ) of the orthorhombic LMO (Fig. 6e). Given a monoclinic angle β PC with four-fold twins (denoted by a ± and b ± , respectively) imposed by the cubic symmetry of the STO substrate (Fig. 6f), the RSMs around the (103) PC reflections result in three split peaks, corresponding to the black crosses in Fig. 6c. Two of them, corresponding to (103) PC and (103) PC reflections, are shifted toward the -L or +L direction at an identical distance in the reciprocal space from the expected pseudotetragonal position estimated from the (00 L) PC reflections, which is related to monoclinic angle β PC . Additionally, the peak at the pseudotetragonal position corresponds to (013) PC and (013) PC reflections. The in-plane reciprocal positions give us information about the inplane a PC -and b PC -axis lattice parameters.
Regarding the unit cell orientation, we further considered the mosaic rotation (δ) of the monoclinic unit cell with respect to the diagonal ([110] C ) direction, as shown in Fig. 6e, to interpret the remaining peaks (red crosses in Fig. 6c). The RSM in Fig. 6c reflects the mosaic tilts leading to the rotation of the peak positions, and further peak splits tracing along the polycrystalline arcs (the yellow lines in Fig. 6c). Because of the mosaic rotation, a clockwise To ensure the model unit cell, we checked the RSM around the (223) PC reflection. According to our models, two variant peaks should appear along the L-axis due to the monoclinic tilt, and two more peaks are created by mosaic rotations with different rotation directions along the polycrystalline arcs ( Supplementary Fig. 4d). Figure 6d displays the resultant reflections around the (223) PC reflection. The black crosses stem from the monoclinic unit cells with no mosaic rotation, and the red crosses correspond to the unit cells subjected to mosaic rotation. It is worth noting that the reciprocal positions of the (223) PC and (223) PC reflections are aligned to the pseudotetragonal L position. This arises from the fact that the two mosaic domains have the mirror symmetry of each other with respect to a {110} PC plane, and should be continuously connected to minimize the elastic energy cost within the domain walls (see the schematics in Fig. 6g). In other words, the originally tilted c-axis vector is rotated so that the in-plane component is parallel to the axis of rotation for matching at the wall. Within the framework of our monoclinic models, we were able to determine the pseudocubic lattice parameters of the LMO film to be a PC = 3.97(1) Å, b PC = 3.84(1) Å, c PC = 3.98(6) Å, β PC = 88.2(1)°, and the mosaic rotation angles to be δ = 0.0(3)°for the upright unit cell and δ = 1.1(7)°for the unit cell subjected to the severe mosaic rotation.

Two types of domain and wall textures
Next, we link the structural information to the surface topography for the two monoclinic unit cells with and without mosaic rotation. To do that, we first try to interpret the orientations of domain walls based on the compatibility relation of ferroelastic domains 50 , Here, we consider two spontaneous strain tensors ε $ ðS 1 Þ and ε $ ðS 2 Þ where strain states S 1 and S 2 come from an identical prototype in high symmetry. And,x is the position vector of the domain wall. According to the quadratic Eq. (9), there are two permissible domain walls between two distinct strain states (e.g., a + and b − ; see Fig. 6g-j, Supplementary Fig. 5 and Note); one is the vertical mirror plane along [110] C called W wall and the other is the inclined plane along [110] C called W' wall. These domain wall orientations in M C -type pseudocubic unit cells are contrasted with domain and wall configurations along <100> PC orientations in rhombohedral PbZr 1-x Ti x O 3 and manganites 51-53 because the M A -type monoclinic pseudocubic unit cells (of which the c-axis tilts toward <110> PC direction) result in domain walls along <100> PC by azimuthal rotation of 45°relative to the present M C -type case. The W wall is formed when the monoclinic tilt directions collide at the wall. To minimize the interfacial elastic energy at the W wall, a + and b − domains are subjected to a clockwise or counterclockwise mosaic rotation (δ = ±1.1(7)°) around the [110] C axis in a symmetrical manner (Fig. 6g, h). For a perfect match at the wall, the mosaic rotation angle should be the monoclinic angle (90 − β PC ) divided by ffiffi ffi 2 p , that is~1.2°. However, this mosaic rotation inevitably makes the interfacial energy between the film and the substrate more unstable. So there is a trade-off problem, and the trade-off will depend on the film thickness. The thickness of the LMO film studied here is~150 nm, but we cannot observe the protruding surface topographic feature when the film is as thin as 35 nm (Supplementary Fig. 6).
Meanwhile, the W' wall is formed when the monoclinic tilt direction continuously changes across the domain wall into a head-to-tail type, as shown in Fig. 6i, j. It doesn't require a mosaic rotation (δ = 0.0(3)°) to relieve the interfacial energy. These expected properties of ferroelastic domain walls show a good agreement with the observed topographic features. We conclude that the regions with the bumpy surface in the AFM images consist of a ± /b ± domains with W walls, reflecting the angle of the slope of the surface, ±1.(0)°( Supplementary Fig. 3), similar to the value of mosaic rotation measured by x-ray diffraction. Meanwhile, the flat regions have W' walls ( Fig. 6j). It is noteworthy that the W' wall is inclined to the substrate by the amount of 56°c onsidering the strain state and the compatibility relation, and the W wall is perpendicular to the substrate (Fig. 6h, j and Supplementary Fig. 5).
The periodicities of stripe domains obtained by AFM and Raman spectromicroscopy are quite different, and they are definitely not one-to-one correspondence (Fig. 1e, g). One is topographically flat and the other has an uneven surface. The flat area (represented by the solid circle in Fig. 1) shows the periodic stripe domains on the Raman map in B 2g mode, while the topographically modulated area (denoted by the dashed circle in Fig. 1) has no distinctive stripes. We estimated that the domain widths of the surface modulated region to be~150 nm (300 nm in periodicity) based on AFM line profile data shown in Supplementary Fig. 3, while the domain widths of the flat region to be~250 nm (~500 nm in periodicity) based on the line profile data from the Raman map shown in Fig. 5b. The seeming inconsistency between the surface topography and the Raman map can be understood based on the spatial resolution of our Raman imaging. The two possible domain textures with different mosaic orientations have the different typical widths of domains. One happens to be larger and the other is smaller than the spatial resolution of our Raman imaging. According to the Kittel's law, the typical size of domains is determined by the competition between the domain energy and the domain wall energy, and the interfacial energy between substrate and film should be also considered. Although it is natural that the typical domain widths for two domain/wall textures are different, further study is required for quantitative theoretical estimates.

Redistribution of the JT domains by mechanical writing
Based on the ferroelastically-formed domain structure, we demonstrate that the distribution of the JT domains could be manipulated by applying a mechanical force, using a AFM tip. We exerted a strong frictional force associated with a normal force of 150 µN corresponding to~100 GPa (approximating the tipsurface contact to be a disk with 20 nm in radius). We scanned along two different directions, [100] C and [110] C , to test whether the trailing effect could align the domains and walls (Fig. 7a). Figure 7 shows the effect of mechanical writing on the topography and Raman signals of the B 2g mode. The mechanical scans erased the W'-wall-type domains (δ = 0°), but the regions became filled with W-wall-type domains (δ = 1.1°). It is attributed to the strong frictional force of the tip, which induces the mosaic rotation and/or forces the monoclinic tilt direction. We note that the written domain walls become denser after the writing. The scan direction also affects the density; the diagonal scan induces a higher density of domain walls with a lower height modulation. The directional control of domain walls doesn't seem to be fully deterministic. The nucleated walls tend to propagate in a straight line, but some have to break the continuation due to the intervention of adjacent configurations. Similar mechanical writing has been accomplished in ferroelastic WO 3 films 48 , and the mechanical load used in the current study is much larger than before. We cannot exclude the possibility that the applied stress would only affect the top surface of the film and thus the depth of the writing would be shallower than the film thickness. Since the scan along [100] C is parallel to the monoclinic tilt direction, the writing effect could reach deeper in the film and the domain width might be wider than the [110] C case.
It is interesting to note that a distinctively enhanced Raman signal arises along the boundaries where the bundles of the Wwall-type domains collide with another orthogonal bundles (indicated by the arrows in Fig. 7c, e). These supermacro domain boundaries are expected to undergo severe local lattice deformation with a large strain and strain gradient and might also be involved in orbital frustration. We note that a strain gradient breaks local inversion symmetry, thereby flexoelectric polarization can be induced 48 . It has been reported that Raman signals are enhanced at ferroelectric domain walls 54 as a consequence of local strain and/or electric fields that can improve Raman scattering efficiency.
We successfully observed the real-space distribution of JT domains in an epitaxial LMO thin film using confocal Raman spectromicroscopy. Based on the symmetry analysis of the Raman modes, we could identify and visualize the local populations and orientation of the JT domains, and the concomitant OO distribution was identified. A detailed structural analysis revealed the formation of ferroelastic domains, which becomes a framework for the JT domain walls under the consideration of the ferroelastic compatibility relation. We also performed mechanical manipulation of the JT domains. Our findings offer a promising approach for the observation of orbital domains, and provide useful information for understanding JT-mediated physics.

Synthesis of the LMO thin films
LMO thin films were deposited on a (001) C -oriented Nb-doped (0.5 wt%) STO substrate (CrysTec GmbH) using pulsed laser deposition with a KrF excimer laser (λ = 248 nm). Growth was made at a temperature of~900°C in the Ar and O 2 gas mixture of 0.1 Torr with a ratio of 9:1. Laser fluence was~0.5 J/cm 2 with a repetition rate of 10 Hz. After deposition, the films were cooled down to room temperature at a rate of 10°C min −1 in vacuum environments (~10 −6 Torr) to resolve the cation vacancy issue, by reducing excess oxygen.

X-ray diffraction analysis of the LMO thin films
To investigate the crystal structure of the LMO film on Nb-doped STO, x-ray θ-2θ scans and RSMs were measured at room temperature using a highresolution x-ray diffractometer (PANalytical X'pert-PRO MRD) with Cu Kα 1 radiation (λ = 1.5406 Å).

Characterization of surface morphology and mechanical writing
The surface topographies of the LMO thin film were characterized using AFM (

Raman spectroscopy
We examined the phonon modes of the LMO thin films using the confocal Raman spectroscope alpha 300R (WITec GmbH, Germany). Raman spectra were excited using a Nd:YAG 532 nm laser, and were recorded in backscattering geometry using a UHTS 600 spectrometer. A laser power density of~20 kW/mm 2 at the sample surface was used. We used an integration time of 0.25 s at each position to obtain the Raman maps of the fundamental phonon modes. To get line profiles of the high-energy regime, an integration time of 30 s was required at each position.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.