Nematicity in the pseudogap state of cuprate superconductors revealed by angle-resolved photoemission spectroscopy

The nature of the pseudogap and its relationship with superconductivity are one of the central issues of cuprate superconductors. Recently, a possible scenario has been proposed that the pseudogap state is a distinct phase characterized by spontaneous rotational symmetry breaking called"nematicity"based on transport and magnetic susceptibility measurements, where the symmetry breaking was observed below the pseudogap temperature $T^*$. Here, we report a temperature-dependent ARPES study of nematicity in slightly overdoped Bi${}_{1.7}$Pb${}_{0.5}$Sr${}_{1.9}$CaCu${}_2$O${}_{8+\delta}$ triggered by a uniaxial strain applied along one of the Cu-O bond directions. While the nematicity was enhanced in the pseudogap state as in the previous studies, it was suppressed in the superconducting state. These results indicate that the pseudogap state is characterized by spontaneous rotational symmetry breaking and that the nematicity may compete with superconductivity. These new experimental insights may provide clues for the nature of the pseudogap and its relation to the superconductivity.


I. INTRODUCTION
The pseudogap in cuprate superconductors is characterized by the suppression of the density of states around the Fermi level below a characteristic temperature T * and above the superconducting transition temperature T c . Broadly speaking, two possible scenarios for the pseudogap have been discussed, that is, a precursor to the superconducting state and distinct order which competes with superconductivity [1]. The latter scenario has been put forward by several measurements which detect distinct orders with spontaneous symmetry breaking such as translational and time-reversal symmetry breaking at or below the pseudogap temperature T * [2][3][4][5][6][7][8][9][10][11][12][13][14][15].
Recently, electronic nematicity that the electronic structure preserves the translational symmetry but breaks the rotational symmetry of the underlying crystal lattice has been found to exist in the pseudogap state [11][12][13][14][15][16][17]. As a possible mechanism, it has been proposed that the nematicity arises from fluctuations of stripe order [18,19] or from the instability of the Fermi surface (so-called Pomeranchuk instability) [20][21][22][23][24]. From experimental perspectives, the nematicity in the cuprate superconductors was first pointed out by transport measurements of lightly-doped La 2-x Sr x CuO 4 and YBa 2 Cu 3 O y (YBCO) [11]. Anisotropic signals in the spin excitation measured by neutron scattering was observed for untwinned underdoped YBCO [12]. Inequivalent electronic states associated with the oxygen atoms in the a and b directions was detected by scanning tunneling spectroscopy measurements of underdoped Bi 2 Sr 2 CaCu 2 O 8+δ (Bi2212) [13]. Nematic fluctuations in Bi2212 have been observed by Raman scattering [25]. Nernst effect and magnetic torque measurements on underdoped and optimally doped YBCO showed a systematic temperature dependence of the nematicity [14,15]. The onset temperature of the in-plane anisotropic signals in the Nernst coefficient and that in the magnetic susceptibility coincide with T * . Furthermore the nematic susceptibility derived from elastoresistance experiments on Bi2212 diverges at T * [26]. The order parameter-like behaviors of the nematicity in the Nernst and magnetic torque experiments [14,15] and the divergence of the nematic susceptibility [26] indicate that the pseudogap state can be considered as a distinct thermodynamic phase characterized by the rotational symmetry breaking. Although the orthorhombic distortion of the CuO 2 plane caused by the Cu-O chains along the b-axis of the untwinned YBCO samples already breaks the four-fold rotational symmetry even above T * , the orthorhombicity is considered to help to form a microscopic nematic domain in one particular direction and enables us to detect nematicity in macroscopic measurements. This weak orthorhombicity of the CuO 2 plane for nematicity plays the same role as a weak external magnetic field for ferromagnets [27].
Motivated by those previous studies on nematicity in the cuprates, we have performed angle-resolved photoemission spectroscopy (ARPES) measurements on slightly overdoped Bi 1.7 Pb 0.5 Sr 1.9 CaCu 2 O 8+δ (Pb-Bi2212) (T c = 91 K) by applying a uniaxial strain along the Cu-O bond direction to detect nematicity below T * . At the doping level we chose, T * was not too high for ARPES experiments but the divergence of the nematic susceptibility in Pb-Bi2212 was sufficiently strong [26]. While most of the studies on nematicity in the cuprate superconductors have been done on YBCO, Pb-Bi2212, whose CuO 2 plane has the tetragonal (C 4 ) symmetry, is a more convenient material to study nematicity than YBCO. Owing to the presence of a natural cleavage plane between the BiO layers and rich information accumulated from previous ARPES studies, Bi2212 is an ideal material for ARPES experiment to investigate novel phenomena in cuprates [28]. In the present work, we analyze the single particle spectra without assuming that the symmetry of the electronic structure is the same as the symmetry of the CuO 2 plane to shed light on the possibly lowered symmetry of the electronic structure. Furthermore, using ARPES, one can investigate the nematicity not only in the normal and pseudogap states but also in the superconducting state, where transport and magnetic measurements cannot be performed due to the zero resistivity and the strong diamagnetism, respectively [14][15][16].

II. EXPERIMENTAL METHODS
Pb-Bi2212 single crystals were grown by the floatingzone method. The hole concentration was slightly overdoped one (T c = 91 K) after annealing the samples in a N 2 flow, which allowed us to measure samples above and below T * (∼160 K) rather easily compared to the underdoped samples whose T * 's are too high [29]. We measured Pb-doped samples in order to suppress the superstructure modulation present in the BiO layers of Bi2212, which causes the so-called diffraction replicas of the Fermi surface shifted by multiples of k = ±(0.21π, 0.21π) [28].
YBCO has Cu-O chains in addition to the CuO 2 plane, which helps nematic domains to align in one particular direction, while Pb-Bi2212 has no such an internal source of strain as the Cu-O chains and the in-plane crystal structure is tetragonal. Therefore, we applied a tensile strain to the sample along the Cu-O bond direction in attempt to align nematic domains in one direction using a device similar to that used for Fe-based superconductors, a picture and a schematic figure of which are shown in Supplementary Fig. S1 (c) [30]. Here, in analogy to ferromagnets, the tensile strain plays a role of a weak external magnetic field that align ferromagnetic domains along the field direction [27]. The strain was applied in air at room temperature, and the stressed sample was introduced into the ultrahigh vacuum and cooled down to measurement temperatures before cleaving in situ (See Supplementary Information S1 for the strain estimated from the lattice parameters by x-ray diffraction measurement). Such an operation was necessary to obtain anisotropic signals from spectroscopic data because the size of the nematic domains is not necessarily larger than the beam size; otherwise, anisotropic signals may be averaged out [31,32]. Note that in our experimental setup, it is impossible to detect nematicity which is diagonal to the crystallographic a and b directions [17], i.e., so-called diagonal nematicity [33] ARPES measurements were carried out at the undulator beamline BL5U of UVSOR using an MBS A-1 analyzer. The photon energy hν was fixed at 60 eV. The energy resolution was set at 30 meV. The linear polarization of the incident light was chosen perpendicular to the analyzer slit and the tilt axis parallel to the analyzer slit, which realizes the unique experimental configuration that preserves the equivalence of the a-and b-axis directions with respect to the light polarization E and the strain direction, as shown in the inset of Fig. 1. This setting guarantees that the anisotropy of the spectroscopic data between the a-and b-axis directions does not originate from matrix-element effects but from intrinsic inequivalence in the electronic structure. The measurements were performed in the normal, pseudogap, and superconducting states in two series, that is, with increasing temperature and with decreasing temperature in order to check the reproductivity. The samples were cleaved in situ under the pressure better than 2 × 10 −8 Pa. Intensity (a. u.) Intensity (a. u.) Intensity (a. u.) Intensity (a. u.) Energy distribution curves (EDCs) of Pb-Bi2212 at various temperatures. The measurements were performed with decreasing temperature for (a)-(d) and increasing temperature for (e)-(g). At each temperature, EDCs at the X and Y points and the difference between them are displayed. Note that the X (Y) point is slightly off from (−π, 0) ((0, π)) to prevent the overlap of diffraction replicas (see also supplementary information S6). To discuss the line shape of the EDCs at each temperature, the spectral weight integrated in the displayed energy range has been normalized. In the pseudogap state (T = 100 K and 125 K), the spectral intensities between the X and Y points are clearly different while they are identical in the normal (T = 200 K and 210 K) and superconducting states (T = 5 K, 25 K, and 50 K). Inset shows the experimental geometry including the crystallographic axes, the strain direction, light polarization, and the analyzer slit, which ensure the spectroscopic equivalence of the a-and baxis directions.

III. RESULTS AND DISUCUSSION
In Fig. 1, energy distribution curves (EDCs) around k = (−π, 0) (X point) and (0, π) (Y point) at various temperatures are shown. The data were obtained from one sample but from two cleavages: on one cleaved surface, we performed ARPES measurements with cooling the sample at T = 5 K, 50 K, 100 K, and 200 K and, on the other cleaved surface, with heating the sample at T = 25 K, 125 K, and 210 K. In both measurements, the tensile strain was applied in the x-direction (||a-axis) as shown in the inset of Fig. 1. In the cooling series, T = 200 K → 100 K → 50 K → 5 K, in the normal state (200 K), the line shapes of the EDCs were almost identical between the X and Y points, as expected from the four-fold rotational symmetry of the CuO 2 plane. With decreasing temperature, the line shapes of the EDCs around the X and Y points became different below T * in the pseudogap state (100 K) and then identical again below T c in the superconducting state (5 K and 50 K). In the heating series, T = 25 K → 125 K → 210 K, the spectral changes of the cooling series was reproduced as shown in Fig. 1. Let us focus on the pseudogap state, where difference in the line shapes of the EDCs is present between the X and Y points. In Bi2212, it is well known that the Cu-O band is split into the anti-bonding and bonding bands due to the bilayer structure [34]. From the dispersions near k = (π, 0) in overdoped sample (T c = 91 K), whose doping level is the same as the present sample, we consider that the energies of the bottoms of the anti-bonding band and bonding band are located around E −E F = -25 meV and -110 meV, respectively, in the pseudogap state. (See Supplementary Information S2.) Thus, we conclude that the intensity of the bonding (anti-bonding) band is higher around the X (Y) point than that around the Y (X) point, reflecting the in-plane anisotropy of the electronic structure in the pseudogap state under the uniaxial strain. Here, we would like to emphasize that this inequivalence between the X and Y points is not caused by matrix-element effect because the x-and y-directions are equivalent for a tetragonal sample in the present measurement geometry (See also Supplementary Information S3).
In Fig. 2 (a)-(g), the intensity maps of the constantenergy surface at E − E F = -30 meV at various temperatures are displayed. We have chosen E − E F = -30 meV rather than E F in the following analysis. This is because the dispersion of quasiparticles near E F is partially gapped due to the pseudogap and superconducting gap opening at low temperatures. The intensity around k y = −k x almost vanishes due to matrix-element effect (See Supplementary Information S4 for the matrixelement effect). Note that additional features which resemble the constant-energy surface but are shifted in k space can be seen in the data of the heating series ( Fig. 2 (e)-(g)). They are due to extrinsic signals aris-ing from a tilted surface that appeared when the sample was cleaved at low temperature. Well-known diffraction replica shifted by (π, π) was also observed. (See Supplementary Information S5 and S6) However, the extrinsic features could be clearly separated from the intrinsic ones in the momentum distribution curves (MDCs), and thus the extrinsic features do not affect the present analysis. At k = (π, 0) for hν = 60 eV, the intensity of the antibonding band is only a little weaker than the bonding band but not negligible [35,36], and their MDC widths are large (∼ 0.3 π/a FWHM) compared to the momentum separation of the two bands (∼ 0.1 π/a). Therefore, we extracted the constant-energy surface practically as a single band rather than multiple bands, i.e., overlapping anti-bonding and bonding bands. The constant-energy surface has been determined from the peak positions of the MDCs fitted using Lorentzians. In Fig. 2 (h), the constant-energy surface at 125 K and 210 K folded in the first quadrant of the first Brillouin zone are overlaid on each other and compared. The constant-energy surface around the Y point at 125 K is closer to (0, π) than that at 210 K. This is qualitatively consistent with the fact that the intensity of the anti-bonding (bonding) band becomes strong near the Y (X) point in the pseudogap state as displayed in Fig. 1 (f). Therefore, we interpret the intensity transfer from the bonding band to anti-bonding band as one goes from the X point to the Y point shown in Fig. 1 (f) makes the constant-energy surface distorted around the Y point.
To understand the anistropy of the constant-energy surface more quantitatively, we have estimated it with the tight-binding model in the following way. The band dispersion near the Fermi level of the cuprates can be fitted using the standard tight-binding model where t, t , and t are the nearest-neighbor, secondnearest-neighbor, and third-nearest-neighbor hopping parameters, respectively, and µ is the chemical potential [37]. The model has four-fold rotational (C 4 ) symmetry.
In order to examine the possibility of the C 4 rotational symmetry breaking into C 2 symmetry, we introduce an anisotropy parameter δ which represents the orthrombicity of the hopping parameters (t and t ) [21] as By fitting the constant-energy surface using Eq. (2), we have estimated its deviation from the C 4 symmetry through the finite δ values (For detailed fitting procedures, see Supplementary Information S7). The temperature evolution of the anisotropy parameter δ thus derived is shown in Fig. 3. In the normal state (T > T * ), δ is close to zero, in both heating and cooling series, which was expected from the four-fold rotational symmetry of the CuO 2 plane. It is also consistent with the EDCs in Figs. 1 (a) and (e). In the pseudogap state (T c < T < T * ), however, δ became finite. The sign of δ indicates that the hopping parameters along the tensile strain direction became small in the pseudogap state. The sign of δ is also the same as the anisotropic Fermi surface of YBa 2 Cu 4 O 8 under the uniaxial strain from the Cu-O chains leading to the small difference between the a and b lattice constants [38]. Thus, the uniaxial strain seems to serve as an external perturbation to align a majority of nematic domains in one direction (See also Supplementary Information S8 for the effect of strain on the magnitude of the hopping parameters). Figure 3 also shows that δ is suppressed in the superconducting state. In the normal and pseudogap states, our result is consistent with the previous magnetic and transport measurements on YBCO in that the nematicity becomes finite below T * [14,15]. More concretely, for example, the anisotropy of the magnetic susceptibility of YBCO, which was derived from the magnetic torque measurements [14] is as large as 0.5% just above T c . Although further studies are necessary to compare the magnetic susceptibility and the single-particle spectral function quantitatively, we naively believe that the δ whose order of the magnitude was 1% in our measurements would have high impacts on various physical properties such as charge and/or magnetic susceptibilities. As for the superconducting state, the magnetic and transport measurements cannot give any information, because of the giant diamagnetism and zero resistivity, respectively, about the underlying electronic structures while our ARPES result indicates possible competition between nematicity and superconductivity.
From the theoretical side, dynamical mean-field theory (DMFT) combined with the fluctuation exchange (FLEX) approximation for the Hubbard model has shown that a Pomeranchuk instability, where the C 2 anisotropy of the Fermi surface is induced, appears in the overdoped region but coexists with superconductivity [37]. In contrast, according to a cellular dynamical mean-field theory (CDMFT) study for the Hubbard model, the C 4 symmetry breaking was shown in the underdoped pseudogap regime rather than in the overdoped regime [39,40]. According to mean-field calculations of the t − J model, the Pomerachuk instability competes with superconductivity [20,41]. Thus, it is not clear at present whether the Pomeranchuk instability consistently explains our result or not.
Nematicity has also been considered as arising from fluctuations of charge-density wave (CDW) or stripes [42,43]. In YBCO, CDW was observed by x-ray scattering experiments and found to reside inside the pseudogap phase and competes with superconductivity [2][3][4][5]. CDW was also found in Bi2212 and to compete with superconductivity [44,45] as in the case of YBCO. More recent work on underdoped Bi2212 (T c = 40 K) showed that the elastic peak that represents the CDW survived up to T * [46]. From the fact that static CDW in Bi2212 competes with superconductivity and possibly survives up to T * including our case, it is also likely that the nematicity observed in our measurements arises from fluctuating CDW. Furthermore, there have been recent theoretical attempts to explain the nematicity using by the pair-density-wave (PDW) model [47,48]. As shown in the calculated single-particle spectra in [48], the inequivalence between the X and Y points may be attributed of the intensity difference predicted for incommensurate PDW states at finite temperatures.
In Fe-based superconductors, nematicity has been observed as the different populations of different orbitals (such as yz vs zx orbitals) revealed by ARPES experiments [31]. In the present study, the observed anisotropy might be due to the different populations of the oxygen p x and p y orbitals because the entire Fermi surface consists of the copper 3d x 2 −y 2 orbital hybridized with the oxygen p x and p y orbitals [49]. However, nematicity can also be triggered by other mechanisms, e.g., the Pomeranchuk instability in single band systems and, therefore, it is not clear whether the anisotropy of the orbital populations is important to explain the nematicity observed in the cuprates as in the case of the Fe-based superconductors. Although at present it is difficult to identify the microscopic origin of the nematicity, our result pro-vides evidence that the pseudogap state shows nematicity and competes with superconductivity. In contrast, the STM experiments have indicated nematicity even in the superconducting state [50]. In order to reconcile the STM experiments with the apparent competition between the nematicity and the superconductivity implied in the present study, there is the possibility that the nematicity is masked by the d-wave superconducting order parameter, whose amplitude should be identical between the two Cu-O bond directions, in the present ARPES data.

IV. CONCLUSIONS
In conclusion, we have demonstrated the presence of nematicity in the pseudogap state of Pb-Bi2212 by temperature-dependent ARPES experiments as suggested in previous studies. On top of that, possible suppression of the nematicity in the superconducting state was also indicated. However, there are still unsolved issues regarding the pseudogap, that is, how the nematicity is related to the other proposed orders inside the pseudogap phase, e.g., CDW [2][3][4][5][6] and loop current order [7][8][9][10]51] which correspond to translational and timereversal symmetry breaking, respectively, and how the Q = 0 nematic order opens pseudogap in the quasiparticle band dispersion. Therefore, further work is necessary to identify the nature and the origin of nematicity. After the ARPES measurements, we carried out x-ray diffraction measurements in order to estimate the magnitude of the uniaxial tensile strain applied to the sample by the screw shown in Fig. S1(c). We estimated the strain from the anisotropy of the a and b lattice parameters because it is difficult to measure the strain directly because of the plasticity of the device. As shown in Fig. S1(a), we focused on the (200) and (020) Bragg peaks of Pb-Bi2212 and each peak was fitted using a single Gaussian. Figure S1 is much smaller than the anisotropy parameter δ of order 1% which we deduced from the ARPES measurements shown in Fig. 3. Therefore, the enhancement of the nematicity in the pseudogap state cannot be explained due merely to the anisotropy of the crystal structure.
A picture and a schematic figures of the strain device are shown in Fig S1(c).

S3. Energy distribution curves in the nordal region
In the nematic state where the electronic structure becomes inequivalent between the k x and k y directions, one expects that the electronic structure around the node is less affected by the nematicity. In order to see if this is the case or not, we compare energy distribution curves (EDCs) around the nodes in the first and third quadrants of the first Brillouin zone at 125 K, the temperature at which the greatest anisotropy was observed between the a and b directions. As shown in Fig. S3, the EDCs around the nodes are almost identical at the two different points in the first Brillouin zone in contrast to the EDCs at the X and Y points shown in the Fig. 1 (F).

S4. Matrix-element effect
In vacuum, the vector potential A( r, t) satisfies the wave equation, and one of the solutions of this equation is where ε denotes the polarization (the oscillating direction of the electric field), and k and ω are the wave vector and the frequency of light, respectively (ω = ck). The probability of photoemission is given by according to Fermi's golden rule. In ARPES measurements, we often use vacuum-ultra violet, whose photon energy is about 100 eV i.e., the wave length is about 100Å which is much longer than the distance between neighboring atoms in the material, so that one obtains e i k· r = 1 + i k · r + · · · 1. (S4.5) Then the matrix element in Eq (S4.3) can be expressed as using the Hermite property of H 0 and the commutation relationship (S4. 10) In summary, the transition probability can be approximated by and the matrix element can be described by that of the dipole transition.
The k-dependence of the intensity in k-space obtained by ARPES can be estimated by calculating the matrix element (Eq. (S4.11)) in the case that |i = |d x 2 −y 2 , E ∝ (+1, −1, 0), which represents our experimental condition, and |f = e i k· r was assumed for simplicity.
Thus, we obtain (S4.14) Figure S4 is the result of the k-dependence of the transition probability in the first Brillouin zone.     Supplementary Figure S11. The same as Fig. S9 except that the data were taken at 25 K.

S6. Diffraction replicas
As seen in the raw MDC data, we have observed some additional features in the constantenergy-surface mapping on top of the intrinsic feature analyzed in the main text. In the cooling series (200 K, 100 K, 50 K, and 5 K), so-called shadow bands, replicas of the main Fermi surface shifted by k = (π, π) of the structural origin [2], are seen. On top of that, the replicas due to a tilted part of the cleaved surface are also seen in the heating series (25K, 125 K, and 210 K). The shift is estimated to be k = (0.08π, 0.25π). Some of those replicas, particularly of the (π, π) replicas, intersect with the main feature in the antinodarl region, where MDC peak positions are crucial to study the distortion of the constant-energy surface representing the nematicity. Therefore, when the intensities of such extrinsic features are not negligibly small, the peak positions of MDCs may be deviated from the intrinsic ones. We have excluded such data points for analyzing the anisotropy of the constant-energy surfaces to eliminate the effect of extrinsic features.

S7. Fitting procedure using the tight-binding model
Fitting procedure to estimate the anisotropy of the constant energy surface shown in Fig.   3 is as follows. First, the experimental data points which represent the constant energy surfaces at 200 K and 210 K, the highest temperature in the cooling and heating series, respectively, were fitted using the tight binding model Eq. (2). Since the actual sample angle, i.e., the small misalignment of the sample, with respect to the sample stage cannot be determined with sufficient accuracy, we chose the the misalignment angles between the sample and the sample stage θ 0 , φ 0 , ω 0 for the rotation around the vertical axis, the tilt, and the azimuth rotation (see Fig. S13), respectively, such that the difference between the data points and fitted curve was minimized. At the highest temperatures in both measurement series (200 K and 210 K), only t /t was fixed with the constraint: −0.51 < t /t < −0.5.
Note that, despite the fact that we did not assume δ = 0 in the procedure to define the sample angles, the anisotropy parameter δ converged to almost zero. At lower temperatures, t /t was also fixed within ± 1 % of t /t at the highest temperatures. The misalignment angles at lower temperatures were defined with the same method as the one at the highest temperatures, however, using Eq. (1), i.e., imposing δ = 0, instead of Eq. (2) where non-zero δ is allowed.
The aim of this further constraint is to avoid the sample angles and the parameter δ from being trapped at a local minimum of the fitting error. After defining the sample angles at each temperature, the data points were fitted with Eq. (2), i.e., δ was allowed to vary.

S8. Effect of strain on hopping parameters
One may suspect that the anisotropy arises only from the effect of the uniaxial strain and that the temperature dependence of the anisotropy is caused by thermal expansion. Such a possibility can be excluded by the following consideration: here, we consider a transfer integral between atomic orbitals of different atoms i and j. When the axis connecting both atoms is taken as the z-axis, the transfer integral between the atoms i and j −t ij (0, 0, 1) = i|h|j is given by the Slater-Koster parameter (l i l j µ), where l A is the azimuthal quantum number of atom A (A = a, b) and µ is the magnetic quantum number of atoms i and j.

respectively.
As for the transfer integral between the 2p orbital of O atom and the 3d orbital of If the anisotropy of t which we observed is only due to the uniaxial strain, the anisotropy parameter δ(= ∆t/t) should change monotonically because the atomic distance changes monotonically with temperature by thermal expansion. However, the anisotropy parameter δ observed in our measurement was enhanced only in the pseudogap state. This means that the stress works to align the nematic domain properly and ensures that the anisotropy does not come only from the distortion of the sample.