## Abstract

Coupling of charge and spin degrees of freedom is a critical feature of correlated electron oxides, as represented by the spin-related mechanism of a Cooper pair under high-*T*_{c} superconductivity. A doublon-holon pair generated on an antiferromagnetic spin background is also predicted to attract each other via the spin-spin interaction *J*, similar to a Cooper pair, while its evidence is difficult to obtain experimentally. Here, we investigate such an excitonic effect by electroreflectance spectroscopy using terahertz electric field pulses in undoped cuprates: Nd_{2}CuO_{4}, Sr_{2}CuO_{2}Cl_{2}, and La_{2}CuO_{4}. Analyses of the spectral changes of reflectivity under electric fields reveal that the splitting of odd-parity and even-parity excitons, a measure of doublon-holon binding energy, increases with *J*. This trend is reproduced by *t-J–*type model calculations, providing strong evidence of the spin-related doublon-holon pairing. Agreement with the calculations supports the s-wave symmetry of the doublon-holon pair in contrast to the d-wave Cooper pair in doped cuprates.

## INTRODUCTION

The coupling of charge and spin degrees of freedom (hereinafter called charge-spin coupling) sometimes creates characteristic physical properties in correlated electron materials. The colossal magnetoresistance caused by the double-exchange mechanism in perovskite manganites and the superconductivity originating from spin fluctuations in high-*T*_{c} cuprates are typical examples (*1*). The charge-spin coupling is also considered to substantially affect the optical properties (*1*–*3*). The electronic structure of undoped cuprates is sometimes discussed using a single-band Hubbard model (*4*). The previous theoretical studies based on a single-band Hubbard model predicted that an electron and a hole or, equivalently, a doublon and a holon photogenerated on a background of antiferromagnetically interacting spins in a Mott insulator attract each other not only via the Coulomb interaction but also via charge-spin coupling, leading to the formation of an excitonic bound state (*5*, *6*). In addition, femtosecond pump-probe spectroscopy revealed that photoexcited doublon-holon pairs decay on the subpicosecond time scale (*7*), and a rapid energy transfer to spin excitations or magnons is suggested to be a plausible mechanism for such an ultrafast decay dynamics of photogenerated charge carriers (*8*–*14*).

Figure 1 shows the schematics of the theoretically predicted spin-related excitonic effect in two-dimensional (2D) Mott insulators (*5*, *6*); the ground state and an excited doublon-holon pair are illustrated in Fig. 1A and Fig. 1, B and C, respectively. In the case where the doublon and holon are apart from each other (Fig. 1B), the energy of a doublon-holon pair is expressed as ∆ + 8*J*. Here, ∆ is the energy of the charge excitation, and *J* is the antiferromagnetic exchange interaction. We define 8*J* as the increase of the energy in the spin sector from the ground state. In the case where a doublon and a holon are located at the neighboring sites, the increase of the energy is ∆ + 7*J* (Fig. 1C). Thus, there should be a binding energy in the order of *J* between a doublon and a holon due to the antiferromagnetic exchange interaction. This picture suggests the spin-related excitonic effect, which is very analogous to the attractive interaction of a Cooper pair in doped cuprates. Presently, however, the influence of charge-spin coupling on the excitonic effect in undoped 2D Mott insulators has not been demonstrated experimentally. In this study, we investigate the origin of the excitonic effect in undoped cuprates by applying terahertz (THz) pulse-pump optical reflectivity probe spectroscopy on three typical compounds, Nd_{2}CuO_{4}, La_{2}CuO_{4}, and Sr_{2}CuO_{2}Cl_{2}, having different magnitudes of *J* of 155, 133, and 125 meV, respectively (*15*, *16*).

Figure 2A shows the crystal structures of three compounds (*17*–*19*). In La_{2}CuO_{4}, each CuO_{2} plane is composed of CuO_{6} octahedrons, while no apical oxygens exist in the CuO_{2} planes of Nd_{2}CuO_{4}. In Sr_{2}CuO_{2}Cl_{2}, the apical oxygens of the CuO_{6} octahedrons are replaced with chlorines. The 2D electron systems of these compounds consist of O 2*p* and Cu 3*d* orbitals. As shown in Fig. 2B, the Mott-Hubbard gap opens in the 3*d* band of Cu because of the large on-site Coulomb repulsion energy *U*. The O 2*p* band is located between the Cu 3*d* upper-Hubbard (UH) band and the lower-Hubbard (LH) band. Thus, the optical gap corresponds to the charge transfer (CT) transition from the O 2*p* valence band to the Cu 3*d* UH band. Figure 2C presents the imaginary part of the dielectric constant, ε_{2}, spectra with the electric fields of lights *E* parallel to the *a* axis (*E*//*a*) in the three compounds (see Results). The energies (*E*_{CT}) of the CT transitions (the peak energies of the ε_{2} spectra) are ~1.6, ~2.2, and ~1.9 eV in Nd_{2}CuO_{4}, La_{2}CuO_{4}, and Sr_{2}CuO_{2}Cl_{2}, respectively.

Electroreflectance (ER), which measures the electric field–induced change of the optical reflectivity, is one of the most effective method to investigate these excitonic effects not only in usual semiconductors (*20*–*22*) but also in correlated electron materials (*23*, *24*). Using this method, we can determine the energies of one-photon–forbidden even-parity states, which cannot be detected by linear optical spectroscopy. In low-dimensional semiconductors such as π-conjugated polymers and silicon polymers, the splitting energy *E*_{s} between the lowest excitonic state with odd parity and the second-lowest excitonic state with even parity can sometimes be a measure of the exciton binding energy in the lowest excitonic state (*21*, *22*) because the second-lowest excitonic state is close to the edge of the electron hole continuum. In a typical ER spectroscopy, optical reflectivity changes are measured by placing two electrodes on both sides (edges) of a crystal (a film) and applying a quasi-static voltage between them. However, with this method, it is difficult to apply strong electric fields on undoped cuprates because they are good semiconductors. Note that a strong electric field generates a large current flow, which leads to the sample degradations. This makes it impossible to measure the reflectivity changes induced by the electric field. To overcome this problem, we use a nearly monocyclic THz pulse (*25*) as an external electric field, in which the duration of the electric field is only ~1 ps. Because of the short pulse duration, this method enables us to apply electric fields far stronger than 100 kV/cm with a negligible current flow. Using such a THz pulse, we can measure the electric field–induced reflectivity changes, which are equivalent to those obtained by the conventional ER spectroscopy.

We report the results of the ER spectroscopy measurement on 2D Mott insulators of three undoped cuprates, Nd_{2}CuO_{4}, La_{2}CuO_{4}, and Sr_{2}CuO_{2}Cl_{2}, using a THz pulse as an electric field source. We can consider an electric field–induced reflectivity change as a kind of third-order nonlinear optical responses and derive the spectrum of the third-order nonlinear susceptibility χ^{(3)} from the results. The analyses of the χ^{(3)} spectra reveal that the even-parity exciton is located at the energy lower than the odd-parity exciton and that the energy difference between them, which is a measure of the exciton binding energy of the former, increases with an increase of *J*. This trend is reproduced by the *t-J*–type model calculations, which provides an unambiguous evidence of spin-related excitonic effects. The agreement of the calculations with the experimental results supports the s-wave symmetry of the excitonic doublon-holon pair in contrast to the d-wave symmetry of a Cooper pair in doped cuprates.

## RESULTS

### Transient reflectivity changes by a THz electric field in Nd_{2}CuO_{4}

First, we show the results in Nd_{2}CuO_{4}. Figure 3A shows the schematic setup of the THz pulse-pump optical reflectivity probe measurement. The polarizations of the THz pump pulse and the optical probe pulse are parallel, and both of them are along the *a* axis. Figure 3B shows a typical electric field waveform, *E*_{THz}(*t*_{d}), of the THz pump pulses. The maximum of ∣*E*_{THz}(*t*_{d})∣ is ~400 kV/cm. The Fourier power spectrum of this pulse is shown by the red line in Fig. 3C. It has a central frequency of 0.8 THz and ranges up to 2.7 THz. The polarized absorption spectrum of a Nd_{2}CuO_{4} single crystal in the THz region with *E*//*a* is also shown by the blue line in the same figure. An absorption peak at ~4 THz is ascribed to an infrared-active phonon mode associated with the displacements of Cu and O ions in the CuO_{2} plane against Nd and O ions in the Nd_{2}O_{2} plane (*26*, *27*). In the lower-frequency region below this absorption, no absorption exists. This result is consistent with the polarized reflectivity spectrum in the THz region previously reported (*26*). Thus, we can consider that the THz pulse does not excite any phonons or electronic transitions in Nd_{2}CuO_{4}.

The reflectivity change Δ*R*(*t*_{d})/*R* at 1.55 eV near the CT transition peak is shown by the open circles in the middle part of Fig. 3D. The reflectivity sharply decreases by the THz electric field around the time origin. The time characteristic of Δ*R*(*t*_{d})/*R* shows exactly the same as that of –[*E*_{THz}(*t*_{d})]^{2} (the blue solid line). This relation Δ*R*(*t*_{d})/*R* ∝ [*E*_{THz}(*t*_{d})]^{2} indicates that the observed response does not originate from a Zener tunneling (*28*) or a Franz-Keldysh effect (*29*, *30*) but is caused by the third-order optical nonlinearity (*31*). At the lower (0.95 eV) and higher (1.70 eV) energies, Δ*R*(*t*)/*R* shows the similar instantaneous changes following [*E*_{THz}(*t*_{d})]^{2}, although the Δ*R*(*t*)/*R* signals are positive.

In the ER spectroscopy using THz electric fields, the photon energy of the THz pulse (~4 meV) was negligibly small compared with the optical gap. In the case that the intervals of the energy levels in the excited states are larger than ~4 meV, THz electric fields can be regarded as static electric fields (*25*). This condition applies to our case, as shown later. Thus, the observed nonlinear optical responses are expressed as follows

Here, *P*^{(3)}(ω) is the third-order nonlinear polarization. *E*_{THz}(0) and *E*(ω) represent the electric fields of the THz pump pulse and the optical probe pulse, respectively. ε_{0} and χ^{(3)}(−ω; 0, 0, ω) are the permittivity of vacuum and the third-order nonlinear susceptibility, respectively.

To obtain information about energy-level structures of the photoexcited state, we measured the probe energy dependence of Δ*R*(*t*_{d})/*R* and plotted the values, Δ*R*(0)/*R*, at the time origin, which are shown in Fig. 3E (open circles), together with the original reflectivity (*R*) spectrum. A characteristic plus-minus-plus structure is observed. Using the Kramers-Kronig (KK) transformation, we obtained the dielectric constant *R* and Δ*R*(0)/*R* spectra, respectively (see the Supplementary Materials). We also calculated the Imχ^{(3)}(−ω; 0, 0, ω) spectrum from the Δε_{2} spectrum using the relationship Δε_{2} = 3[*E*_{THz}(0)]^{2} Imχ^{(3)}(−ω; 0, 0, ω). The obtained ε_{2} and Imχ^{(3)}(−ω; 0, 0, ω) spectra are shown in Fig. 4A.

The Imχ^{(3)}(−ω; 0, 0, ω) spectrum also shows a plus-minus-plus structure around the CT transition peak. This spectrum cannot be reproduced by a simple narrowing or broadening of the original ε_{2} spectrum (see the Supplementary Materials). Imχ^{(3)}(−ω; 0, 0, ω) or Δε_{2} spectra in 1D semiconductors previously obtained by conventional ER spectroscopy revealed that a one-photon–forbidden even-parity state ∣2〉 was located above a one-photon–allowed odd-parity state ∣1〉 and that the field-induced mixing of those two states caused a plus-minus-plus structure (*20*–*22*). In this case, the ε_{2} and Imχ^{(3)}(−ω; 0, 0, ω) or Δε_{2} spectra can be analyzed using the three-level model, in which the ground state ∣0〉 and two excited states ∣1〉 and ∣2〉 are considered. In our study, we also adopt a similar three-level model.

Before the analyses of the Imχ^{(3)}(−ω; 0, 0, ω) spectra using the three-level model, it is necessary to analyze the original ε_{2} spectrum. To express the dispersion of ε_{2} originating from the odd-parity state, we adopt the following Lorentzian-type dielectric function*N* is the density of Cu atoms, *e* is the elementary charge, and ℏ is the reduced Planck constant. 〈0∣*x*∣1〉 represents the transition dipole moment between ∣0〉 and ∣1〉. ℏω_{1} and ℏγ_{1} are the energy and the damping constant of the state ∣1〉, respectively. We evaluated the magnitudes of 〈0∣*x*∣1〉, ℏω_{1}, and ℏγ_{1} by fitting Eq. 2 to the experimental ε_{2} spectrum. The fitting curve (the red dashed line at the top of Fig. 4A) reproduced well the peak structure in the ε_{2} spectrum.

Assuming the three-level model, we next simulate the Imχ^{(3)}(−ω; 0, 0, ω) spectrum. According to the perturbation theory, χ^{(3)}(−ω; 0, 0, ω) consists of 12 terms (*31*). The main term is expressed as

Here, 〈1∣*x*∣2〉 is the transition dipole moment between ∣1〉 and ∣2〉, and ℏω_{2} and ℏγ_{2} are the energy and damping constant of state ∣2〉, respectively. Using the complete expression of χ^{(3)}(−ω; 0, 0, ω) including 12 terms, we simulated the experimental Imχ^{(3)}(−ω; 0, 0, ω) spectrum. In this analysis, we used the values of ℏω_{1}, ℏγ_{1}, and 〈0∣*x*∣1〉 obtained from the analysis of the ε_{2} spectrum. Thus, the fitting parameters were ℏω_{2}, ℏγ_{2}, and 〈1∣*x*∣2〉. The red dashed line at the bottom of Fig. 4A shows the fitted curve, which reproduces the experimental spectrum well. The energy positions of the odd-parity and even-parity states (ℏω_{1} and ℏω_{2}) are indicated by the dotted lines in Fig. 4A. The successful fitting suggests that the states ∣1〉 and ∣2〉 have discrete energy levels and are excitonic bound states or equivalently excitons consisting of a doublon-holon pair. The even-parity exciton is located at the lower energy than the odd-parity exciton.

As seen in Fig. 4A, the higher-energy part of the ε_{2} spectrum colored in green cannot be reproduced by the Lorentzian function (Eq. 2). From the spectral shape of this component, it is reasonable to consider that it corresponds to the continuum and that the odd-parity exciton responsible for the peak of the ε_{2} spectrum is located close to the band edge. Thus, the splitting *E*_{s} = ℏω_{1} − ℏω_{2} between the odd-parity and even-parity excitons can be a measure of the binding energy of a doublon-holon pair in the lowest even-parity exciton. The origin of the excitonic effect is elucidated in Discussion.

The values of the physical parameters evaluated from the analyses of the linear and nonlinear optical spectra, ℏω_{1}, ℏω_{2}, ℏγ_{1}, ℏγ_{2}, 〈0∣*x*∣1〉, and 〈1∣*x*∣2〉, the maximum values of ∣Imχ^{(3)}(−ω; 0, 0, ω)∣[max∣Imχ^{(3)}(−ω; 0, 0, ω)∣], and *E*_{s} are listed in Table 1, together with the other important parameters, the exchange interaction *J* and the Cu─O bond length *d*_{Cu─O} in the 2D CuO plane.

### Transient reflectivity changes by a THz electric field in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}

To investigate how the nonlinear response to the THz electric field observed in Nd_{2}CuO_{4} depends on physical parameters, such as the optical gap energy *E*_{CT} (ℏω_{1}), the exchange interaction *J*, and the structure of the CuO_{2} planes, we performed similar studies in the other cuprate Mott insulators, La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}. The previous polarized reflectivity measurements in the THz region revealed that the lowest infrared-active phonons are located at 145 and 176 cm^{−1} in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}, respectively, and that no absorptions exist below 100 cm^{−1} (*27*), so that the reflectivity changes generated by a THz pulse in these two compounds can be considered to be the nonlinear electronic responses to the THz electric field, similar to the case of Nd_{2}CuO_{4}.

In La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}, we measured the time evolutions of reflectivity changes, Δ*R*(*t*_{d})/*R*, by the THz pulse in almost the same experimental condition as in Nd_{2}CuO_{4}. The maximum electric field of the THz pulse is ~400 kV/cm. The detected Δ*R*(*t*_{d})/*R* signals follow the square of the THz electric field, [*E*_{THz}(*t*_{d})]^{2}. The spectra of Δ*R*(0)/*R* show a plus-minus-plus structure similar to that of Nd_{2}CuO_{4} in Fig. 3E. The data in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2} are reported in the Supplementary Materials.

We performed the KK transformations of the Δ*R*(0)/*R* spectra and obtained the Imχ^{(3)}(−ω; 0, 0, ω) spectra, which are shown in Fig. 4 (B and C), together with the ε_{2} spectra. The Imχ^{(3)}(−ω; 0, 0, ω) spectra of Sr_{2}CuO_{2}Cl_{2} and La_{2}CuO_{4} also show a plus-minus-plus structure. We analyzed the ε_{2} and Imχ^{(3)}(−ω; 0, 0, ω) spectra with the same procedures performed on those spectra of Nd_{2}CuO_{4} mentioned in the previous subsection. The fitting curves (Eq. 2) of the lowest peak structures in the ε_{2} spectra are shown by the red dashed lines at the top of Fig. 4 (B and C). The peaks are attributable to the odd-parity excitons, and the deviations of the data from the fitting curves colored in green are also considered to correspond to the doublon-holon continuum. The Imχ^{(3)}(−ω; 0, 0, ω) spectra are analyzed by the three-level model. The fitting curves are shown by the red dashed lines at the bottom of Fig. 4 (B and C), which almost reproduce the characteristic plus-minus-plus structures of Imχ^{(3)}(−ω; 0, 0, ω). The parameter values evaluated from the analyses of the linear and nonlinear optical spectra in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2} are also listed in Table 1, together with *J* and *d*_{Cu─O} values.

As seen in Fig. 4 (A to C), in all the three compounds, the energy of the even-parity exciton ∣2〉 was lower than that of the odd-parity exciton ∣1〉, i.e., *E*_{s} = ℏω_{1} − ℏω_{2} > 0. The even-parity excitons in 2D cuprate Mott insulators were previously studied by two-photon absorption spectroscopy (*32*, *33*). However, the reported energy positions depend on the experimental conditions. In pump-probe absorption spectroscopy of Sr_{2}CuO_{2}Cl_{2}, the even-parity and odd-parity states are almost degenerate (*32*), whereas in Z-scan spectroscopy of Nd_{2}CuO_{4}, the even-parity state is located below the odd-parity state (*33*). The present study supports the latter result. In the 1D Mott insulators such as halogen-bridged nickel-chain compounds, the even-parity exciton is located at the higher energy than the odd-parity exciton (*E*_{s} = ℏω_{1} − ℏω_{2} < 0) (*23*). This feature is the same as that in the other 1D semiconductors such as π-conjugated polymers and silicon polymers. The feature of 1D Mott insulators is that the difference between ℏω_{2} and ℏω_{1} is very small. In a bromine-bridged nickel-chain compound, [Ni(chxn)_{2}Br]Br_{2} (chxn, cyclohexanediamine), it is just 10 meV, while the optical gap energy is ~1.3 eV (*24*), that is, the odd-parity and even-parity excitons are almost degenerate in 1D Mott insulators. This is related to the nature of spin-charge separation characteristic of 1D correlated electron systems; the charge sector is separated from the spin sector, and the wave functions of the odd-parity and even-parity excitons are almost the same with each other except for their phases (*34*). Therefore, the lower-lying even-parity exciton demonstrated in the present study on the three undoped cuprates is considered to be characteristic of 2D Mott insulators (*35*). In Fig. 5A, we show the energy-level structure of 2D cuprate Mott insulators. The symmetry of the wave functions of two excitonic states is discussed again later.

## DISCUSSION

First, we discuss the material dependence of the parameters obtained from the analyses of the linear and nonlinear optical spectra in the three-cuprate Mott insulators. To examine the trends of the physical parameters, we plot in Fig. 5 (B to D) the values of *J*, *E*_{s}(= ℏω_{1} − ℏω_{2}), max∣Imχ^{(3)}(−ω; 0, 0, ω)∣, 〈0∣*x*∣1〉^{2}, and 〈1∣*x*∣2〉^{2} as a function of the lowest odd-parity exciton energy ℏω_{1}. The parameters except for 〈0∣*x*∣1〉^{2} seem to be closely correlated with each other.

According to the simple perturbation treatment, 〈0∣*x*∣1〉 is proportional to *t*_{pd}/*E*_{CT} (*36*), where *E*_{CT} = ℏω_{1} and *t*_{pd} is the transfer energy between the O 2*p* and the Cu 3*d* orbitals. We can consider that *t*_{pd} is dominated by *d*_{Cu─O}; *t*_{pd} increases with the decrease of *d*_{Cu─O}. The relatively large 〈0∣*x*∣1〉 in Na_{2}CuO_{4} is ascribed to the small ℏω_{1} (Table 1). La_{2}CuO_{4} has the largest ℏω_{1}, but its 〈0∣*x*∣1〉 is slightly larger than that of Sr_{2}CuO_{2}Cl_{2}. This is attributed to the relatively small *d*_{Cu–O} (large *t*_{pd}) in La_{2}CuO_{4}. However, the difference in 〈0∣*x*∣1〉 in the three compounds is not so large (Table 1). The similar perturbation treatment gives a relation *J* ∝(*t*_{pd}/*E*_{CT})^{2} (*36*). The subtle balance of *t*_{pd} and *E*_{CT} leads to the largest *J* (155 meV) in Na_{2}CuO_{4}, the smallest *J* (125 meV) in Sr_{2}CuO_{2}Cl_{2}, and the intermediate *J* (133 meV) in La_{2}CuO_{4}.

Among the three compounds, Sr_{2}CuO_{2}Cl_{2} has the largest max∣Imχ^{(3)}(−ω; 0, 0, ω)∣. The important factors that dominate the magnitude of ∣Imχ^{(3)}(−ω; 0, 0, ω)∣ are the transition dipole moments, 〈0∣*x*∣1〉 and 〈1∣*x*∣2〉, because ∣Im χ^{(3)}(−ω; 0, 0, ω)∣ is proportional to 〈0∣*x*∣1〉^{2}〈1∣*x*∣2〉^{2}. Figure 5D shows that the relative magnitudes of max∣Imχ^{(3)}(−ω; 0, 0, ω)∣ are dominated by those of 〈1∣*x*∣2〉^{2}; the largest max∣Im χ^{(3)}(−ω; 0, 0, ω)∣ in Sr_{2}CuO_{2}Cl_{2} is caused by the largest 〈1∣*x*∣2〉^{2}, which is much larger than that in the other two compounds.

In Fig. 5 (C and D), 〈1∣*x*∣2〉^{2} seems also related to *E*_{s}: 〈1∣*x*∣2〉^{2} increases from Nd_{2}CuO_{4} to Sr_{2}CuO_{2}Cl_{2} and then decreases to La_{2}CuO_{4}, while *E*_{s} decreases and then increases, leading to the largest 〈1∣*x*∣2〉^{2} and the smallest *E*_{s} in Sr_{2}CuO_{2}Cl_{2}. Such an inverse correlation between 〈1∣*x*∣2〉^{2} and *E*_{s} can be understood if one notices the fact that a small *E*_{s} implies similar spatial extensions of the wave functions along the *a* axis in the states ∣1〉 and ∣2〉, which may lead to a large 〈1∣*x*∣2〉. *E*_{s} in Fig. 5C does not depend systematically on the optical gap energy ℏω_{1} but shows a strong correlation with *J* in Fig. 5B. By plotting *E*_{s}, 〈0∣*x*∣1〉^{2}, 〈1∣*x*∣2〉^{2}, and max∣Imχ^{(3)}(−ω, 0, 0, ω)∣ as a function of *J* in Fig. 5 (E and F), we notice that *E*_{s} is scaled by *J* and 〈1∣*x*∣2〉^{2} is also correlated with *J*.

In Results, we pointed out that the magnitude of *E*_{s} can be a measure of the binding energy of the lowest even-parity exciton. The strong correlation between *E*_{s} and *J* makes us expect an important influence of the spin sector on the excitonic effect as we commented on in the introduction referring to Fig. 1. To demonstrate this scenario, we theoretically analyzed the photoexcited states in 2D Mott insulators, where a doublon and a holon are created by an optical pulse and where spins in the background interact with each other by the nearest neighbor exchange interaction *J*. The resulting Hamiltonian for the two particles resembles a *t-J*–type model with two-hole carriers. More precisely, the Hamiltonian contains the exchange interaction *J* and hopping terms of the holon and doublon with the nearest neighbor amplitude *t*, the second-nearest neighbor amplitude *t*′, and the third-nearest neighbor amplitude *t*″, together with the nearest neighbor attractive Coulomb interaction ( −*V*) between the holon and the doublon (*37*). The parameters *t*′ and *t*″ are necessary to introduce the electron hole asymmetry.

The parameters are determined as follows. In the cuprates, it is known that *J*/*t* is almost equal to 0.4 in common (*38*–*40*). Therefore, we set *J*/*t* = 0.4. For *t*′ and *t*″, we use typical relations in the cuprates, *t*′/*t* = −0.25 and *t*″ /*t* = 0.12, respectively (*41*). Because we want to know the influence of the spin sector on the excitonic effect, we use *J* as a parameter to control the sector. To clarify the effect of *V* on the excitonic effect, we investigated *E*_{s} for the three cases, *V*/*t* = 0, 1.5, and 3. By numerically diagonalizing the Hamiltonian for a *6*). Therefore, the eigenenergy of the even-parity state corresponds to ℏω_{2}, whereas that of the odd-parity state corresponds to ℏω_{1}. Figure 5G shows the calculated *E*_{s} as a function of *J* for three values of *V*/*t*. *E*_{s} increases with the increase of *J*, which is qualitatively consistent with the experimental results (Fig. 5E). This suggests that the attractive force between doublons and holons would be related to the energy change of the spin sector. With the increase of *V*, both ℏω_{1} and ℏω_{2} are expected to decrease, while *E*_{s} does not depend on *V* so much (Fig. 5G). This suggests that *J* increases the stability of the even-parity exciton relative to the odd-parity exciton corresponding to an absorption peak that is close to the edge of the doublon-holon continuum. To reproduce the experimental results more quantitatively, an advanced theoretical analysis, e.g., dynamical mean field theory calculation, will be necessary. Such a study is an important subject in the future.

Encouraged by the possible description of the observed excitonic effects using the *t*-*J–*type model with doublon and holon, we discuss the symmetry of the lowest even-parity exciton and the second-lowest odd-parity exciton obtained by the same model. The odd-parity exciton should have a p-wave symmetry (Fig. 5A). The theoretical study also indicates that the lowest even-parity exciton has an s-wave symmetry (Fig. 5A) (*6*). This disagrees with the d-wave symmetry of a Cooper pair in doped cuprates (*42*). A Cooper pair is formed by two holons, whereas an exciton is formed by a doublon and a holon. The doublon has a negative charge, in contrast to the positively charged holon. This charge difference induces an additional sign for the nearest neighbor hopping of a doublon. The sign affects the phase of the pair wave function, leading to s-wave symmetry for a pair of a doublon and a holon in contrast to the d-wave symmetry for two holons.

In summary, we evaluated the χ^{(3)}(−ω, 0, 0, ω) spectra for the 2D Mott insulators of undoped cuprates by the ER spectroscopy using THz electric field pulses. Analysis with a three-level model indicated that the even-parity exciton is located below the odd-parity exciton and that their energy difference, which is a measure of the excitonic effect, increases with increase of the antiferromagnetic exchange interaction *J*. This trend was reproduced by the *t-J*–type model calculations. Our results demonstrate that a doublon and a holon in a 2D Mott insulator attract each other via the antiferromagnetic exchange interaction. As reported here, nonlinear optical spectroscopy using a strong THz pulse is a powerful tool to investigate the hidden nature in electronic structures and photoexcited states. We expect that it can be applied not only to correlated electron systems but also to the other various kinds of complex electron materials.

## MATERIALS AND METHODS

### Sample preparation

In the present study, we used single crystals of Nd_{2}CuO_{4}, La_{2}CuO_{4}, and Sr_{2}CuO_{2}Cl_{2}, which were grown by the methods previously reported in (*43*), (*44*), and (*19*), respectively. The crystal structures of Nd_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2} are tetragonal, and their *a* and *b* axes are parallel to the Cu─O bonds in the CuO_{2} planes. On the other hand, the crystal structure of La_{2}CuO_{4} is orthorhombic in the strict sense (*17*). Here, the axes of La_{2}CuO_{4} are taken in the same manner as Nd_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}, i. e., *a* and *b* axes are parallel to the Cu─O bonds in the CuO_{2} plane as shown in Fig. 2A.

### THz time-domain spectroscopy

To measure the absorption spectrum of Nd_{2}CuO_{4} shown in Fig. 3C, we used a THz time-domain spectroscopy (THz-TDS). The thickness of the sample was 25 μm. As a light source of the THz-TDS, we used a Titanium:sapphire regenerative amplifier (RA) with a pulse duration of 25 fs, a central photon energy of 1.56 eV, a repetition rate of 1 kHz, and an output power of 2.5 W. The output of the RA was divided into two beams. One was used to generate THz pulses via an air plasma–induced THz radiation (*45*). The other was used for the detection of the THz pulses by an electro-optical (EO) sampling. We measured temporal waveforms of the THz pulses transmitted through the sample [*E*_{sample}(*t*_{d})] and the air [*E*_{ref}(*t*_{d})] by the EO sampling with a 300-μm-thick GaP crystal. From the Fourier transformation of the obtained waveforms, the power spectrum of each transmitted THz pulse, *I*_{sample}(ω) and *I*_{ref}(ω), was obtained. Using these spectra, we simply calculated the optical density (OD) spectrum by OD = −log[*I*_{sample}(ω)/*I*_{ref}(ω)].

### THz pulse-pump optical reflectivity probe spectroscopy

As the light source of the pump and probe pulses, we used a Titanium:sapphire RA with a pulse duration of 90 fs, a central photon energy of 1.56 eV, a repetition rate of 1 kHz, and an output power of 4.5 W. The output of the RA was divided into two beams. One was used to generate THz pulses via the optical rectification in a LiNbO_{3} crystal. To generate intense THz pulses, we used a pulse front–tilting method (*46*). The THz pump pulses were focused on the sample by parabolic mirrors. To measure the electric field waveforms, we used a ZnTe crystal instead of a sample and applied an EO sampling method. Detailed procedures for creating and detecting THz pulses are reported elsewhere (*25*). The other output from the RA was used to excite an optical parametric amplifier from which light pulses with 0.4 to 2.8 eV having a temporal width of 90 fs were obtained. The delay time *t*_{d} of the probe pulse relative to the pump pulse was controlled by changing the path length of the probe pulse. The time origin was defined at the position of the maximum THz electric field strength. All the measurements were performed at 294 K. To avoid sample degradation, we performed the measurements of Sr_{2}CuO_{2}Cl_{2} in vacuum.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/6/eaav2187/DC1

Section S1. KK analyses of reflectivity spectra

Section S2. Electric field dependence of electric field–induced reflectivity changes

Section S3. Analyses of electric field–induced reflectivity changes

Section S4. Transient reflectivity changes by the THz pulse and their analyses in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}

Section S5. Contributions of continuum states on Reχ^{(3)}(−ω; 0, 0, ω) and Imχ^{(3)}(−ω; 0, 0, ω) spectra in Sr_{2}CuO_{2}Cl_{2}

Fig. S1. Comparison of polarized absorption (α) spectra with the electric fields of lights parallel to the *a* axis.

Fig. S2. Electric field dependence of reflectivity changes Δ*R*(0)/*R* at the time origin in Sr_{2}CuO_{2}Cl_{2}.

Fig. S3. Reχ^{(3)}(−ω; 0, 0, ω) and Imχ^{(3)}(−ω; 0, 0, ω) spectra and their analyses in Nd_{2}CuO_{4}.

Fig. S4. Transient reflectivity changes in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}.

Fig. S5. Reχ^{(3)}(−ω; 0, 0, ω) and Imχ^{(3)}(−ω; 0, 0, ω) spectra and their analyses in La_{2}CuO_{4} and Sr_{2}CuO_{2}Cl_{2}.

Fig. S6. Comparison of the fitting analyses of Reχ^{(3)}(−ω; 0, 0, ω) and Imχ^{(3)}(−ω; 0, 0, ω) spectra with three-, four-, and five-level models in Sr_{2}CuO_{2}Cl_{2}.

Table S1. Parameter values with errors obtained from the fitting analyses on the ε_{2} spectra.

Table S2. Parameter values with errors obtained from the fitting analyses on the Reχ^{(3)}(−ω; 0, 0, ω) and Imχ^{(3)}(−ω; 0, 0, ω) spectra.

Table S3. Parameter values obtained from the fitting analyses with three-, four-, and five-level models in Sr_{2}CuO_{2}Cl_{2}.

Table S4. Products of dipole moments associated with nonlinear optical spectra in Sr_{2}CuO_{2}Cl_{2}.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was partly supported by a Grant-In-Aid (KAKENHI: project number JP25247049) from the Ministry of Education, Science, Culture, and Sports of Japan and by CREST (grant number JPMJCR1661), Japan Science and Technology Agency. T.Te. was supported by JSPS through the Program for Leading Graduate Schools (MERIT).

**Author contributions:**T.Te., T.O., T.Mi., T.Mo., and H.Y. constructed the THz pump optical probe spectroscopy system. T.Te., T.Mi., and N.K. constructed the THz-TDS system. T.Te. and T.O. performed the measurements. T.I. provided single crystals of Nd

_{2}CuO

_{4}and La

_{2}CuO

_{4}. T.S. provided single crystals of Sr

_{2}CuO

_{2}Cl

_{2}. T.To. performed the theoretical calculations. H.O. coordinated the study. All authors discussed the results and contributed to writing the paper.

**Competing interests:**The authors declare that they have no competing financial interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).