Magnetically brightened dark electron-phonon bound states in a van der Waals antiferromagnet

In van der Waals (vdW) materials, strong coupling between different degrees of freedom can hybridize elementary excitations into bound states with mixed character1–3. Correctly identifying the nature and composition of these bound states is key to understanding their ground state properties and excitation spectra4,5. Here, we use ultrafast spectroscopy to reveal bound states of d-orbitals and phonons in 2D vdW antiferromagnet NiPS3. These bound states manifest themselves through equally spaced phonon replicas in frequency domain. These states are optically dark above the Néel temperature and become accessible with magnetic order. By launching this phonon and spectrally tracking its amplitude, we establish the electronic origin of bound states as localized d–d excitations. Our data directly yield electron-phonon coupling strength which exceeds the highest known value in 2D systems6. These results demonstrate NiPS3 as a platform to study strong interactions between spins, orbitals and lattice, and open pathways to coherent control of 2D magnets.


Report for "Magnetically brightened dark electron-phonon bound states in a van der Waals antiferromagnet".
In this manuscript, Ergecen and Ilyas et al revealed a unique electron-phonon bound state below the Néel transition in NiPS3, and more importantly, identified the electronic origin as localized d-d excitations and the phonon source as A1g phonon at ~ 7THz, using transient absorption spectroscopy with a dynamic range. Van der Waals (vdW) magnets emerge as a flourish platform that hosts strong coupling between multiple degrees of freedom and realizes novel bound states down to the two-dimension (2D) limit. Yet, it is highly challenging to detect the presence of these bound states and more so to resolve the composition of them. This manuscript has shown, for the first time in vdW magnets, the unambiguous identification of electronic and phononic origins of an electron-phonon bound state with the strongest electron-phonon coupling known so far. This manuscript is of high interest and high quality. I would recommend it for publication at Nature Communications, after addressing the following questions.
1. Could the authors please provide an equilibrium absorption spectrum? How does this equilibrium spectrum compare with the transient absorption spectra? Are the fine oscillations for the electronphonon coupled states resolvable in the equilibrium spectrum?
2. What is the role of demagnetization with the high intensity pump in detecting the phonon replicas in the transient absorption spectroscopy? Why is the transient reflectivity probed 2 picoseconds after the pump? What is the temporal dynamics of the phonon replicas?
3. When fitting the spectrum with a sum of Gaussians weighted by a Poisson distribution to extract the Huang-Rhys factor g, there are a few questions that I am interested in: a. Why do the authors choose Gaussian, instead of Lorentzian? b. How do the authors determine where the zeroth order of the Poisson distribution is? c. Is there a broad background in the absorption spectrum? What happens if fitting with a broad background plus a sum of Gaussians weighted by a Poisson distribution? d. The linewidth of individual Gaussians should be the convolution between the phonon linewidth and the d-d excitation linewidth. However, the d-d excitation linewidth, reported in PRL 120, 136402 (2018), seems to be much larger than the Gaussian linewidth here. Any reasons?

Response to Reviewers
The manuscript by Ergecen et al. reported the observation of phonon replicas through the ultrafast optical spectroscopic study of antiferromagnetic vdW insulator NiPS3. These phonon replicas only appeared below Neel temperature. So the authors called these states as magnetically brightened dark electron-phonon bound states. Furthermore, the authors employed the energy resolved coherent phonon spectroscopy to differentiate the origin of these electron-phonon bound states. They found that the coupling between localized d-d transitions and A1g phonon mode is more relevant, in comparison to the coupling between spin-orbit-entanged excitons and A1g phonon mode. The results are interesting and may be considered for the publication in Nature Communications.
Our response: We thank the reviewer for carefully reading our work and providing insightful remarks, and recommending for publication. Here we try to provide comprehensive answers to the questions.
Here I have some technical questions about the work.
1. There are many phonon modes observed in Raman spectroscopy, as seen in ref. 11. And the A1g phonon mode at ~7 THz or ~225 cm-1 is not so prounced in Raman spectra. Why is this specific phonon mode coupled to the d-d transitions?
Our response: First, we would like to clarify the phonon mode that strongly couples to the d-d transitions and gives rise to replica formation. As shown in the Supplementary Information, Fig. S5, the Fourier transform of the replica signal is broad and spans an energy interval between 25 meV to 35 meV. This range includes both the A1g phonon mode with 253 cm-1 wavenumber (31 meV) and Eg phonon mode with 225 cm-1 wavenumber (28 meV). Since our transient absorption measurements do not have enough energy resolution to distinguish between these phonons, we perform energy-resolved coherent phonon spectroscopy, which has a better frequency resolution than our transient absorption measurements, to pinpoint the phonons that couple to the d-d transitions.
In NiPS3, as observed by our coherent phonon spectroscopy measurements, d-d transitions couple to three distinct phonon modes: 1) 5.2 THz phonon oscillation, corresponding to an Eg phonon mode with 173 cm-1 wavenumber. This phonon mode has been reported in Raman measurements (Scientific Reports 6, 20904 (2016)). 2) 7.5 THz phonon oscillation, corresponding to an A1g phonon mode with 253 cm-1 wavenumber. This phonon mode has been denoted as P5 in Ref. 11, and carries a strong spectral weight as evidenced in Raman measurements. In the manuscript, this phonon mode was referred to as "~7 THz A1g mode". 3) 11.5 THz phonon oscillation, corresponding to an A1g phonon mode with 384 cm-1 wavenumber. This phonon mode has been reported in Raman measurements (Scientific Reports 6, 20904 (2016)).
The Eg phonon mode at 225 cm-1 which has negligible spectral weight in Raman measurements does not appear in our coherent phonon spectroscopy as shown in our Supp. Fig. S7, and therefore has negligible coupling to the d-d transitions. In the light of our coherent spectroscopy measurements, the energy of the A1g phonon mode with 253 cm-1 wavenumber (31 meV) matches the spectral distance between the replicas, as shown in the Fourier transform of the phonon replicas (Supp. Fig. S5).
The energy splitting between the d-d levels is dictated by the overlap between the ligand p-orbitals and d-orbitals. Any modulation that alters this overlap (such as strain, pressure and phonon excitation) will cause a shift in d-d energies. As the reviewer pointed out, in our measurements we only observe phonon replicas formed by the d-d transition and A1g phonon mode with 253 cm-1 wavenumber (31 meV). As shown in Figure 4 of Scientific Reports 6, 20904 (2016), this A1g phonon mode corresponds to the out-of-plane motion of sulfur atoms. This distortion modulates the distance between local nickel sites and sulfur ligands, and therefore leads to replica formation by modulating the d-d transition energy. On the other hand, the other A1g phonon mode with 384 cm-1 wavenumber (11.5 THz) corresponds to the collective out-of-plane motion of nickel and sulfur sites with a slight in-plane component for the sulfur sites. Since this phonon mode does not significantly alter the ligand-transition metal distance, it does not contribute to the replica formation.
We thank the reviewer for insightful comments, and we have made the following changes in the text to clarify the phonon modes that couple to the d-d transition and their signatures in Raman and coherent phonon spectroscopy:

A) [They manifest themselves… -Line 38]
-We clarified the frequency of the phonon mode that participates in the replica formation.

E) [We observe two … -Line 94]
-We added a new sentence that summarizes the phonon modes observed in our coherent phonon spectroscopy measurements. F) [The frequency of the dominant … -Line 95] -We clarified the frequency of the phonon mode that participates in the replica formation. G) [These spectral features start … -Line 130] -We clarified the frequency of the phonon mode that participates in the replica formation.
2. The stripy phase of antiferromagnetic order in NiPS3 is globally inversion symmetric. To understand the correlation between phonon replica and magnetic order, the authors proposed a picture of local inversion symmetry breaking. I feel uneasy about this picture. There may be other possibilities such as magnetic point defects and stacking faults in bulk crystals. I suggest the authors add their comments in the manuscript.
Our response: Because of dipole transition rules, transitions between d-levels are not allowed if the inversion symmetry of the transition metal atom is not broken. In the case of NiPS3, the d-d transitions are silent above the magnetic ordering temperature, as the nickel atom is an inversion center. The appearance of d-d transitions at the onset of magnetic order requires a mechanism that links the long range magnetic order and the on-site inversion symmetry breaking. The stacking faults would not give rise to temperature dependent d-d transitions and replica peaks. In addition, magnetic point defects and dislocations can cause loss of inversion symmetry adjacent to the defect sites and can couple to the magnetic order. However, even though defects can influence the optical properties in their vicinity, we think that they cannot give rise to an optical signal that is independent of sample position. In addition, existence of sharp magnetic excitons in our samples is also indicative of high sample quality and sparse defect distribution.
Following the reviewer's comments, we have added the following comments in the revised text, stating that magnetic point defects and stacking faults are unlikely to result in temperature dependent phonon replicas of d-d transitions: 1) [Although stacking faults and lattice defects... -Line 138] -We added a sentence to discuss the stacking faults and lattice defects that can give rise to inversion symmetry breaking.
3. A bulk crystal of NiPS3 cannot be called 2D magnet. Can the authors show the result in few-layer NiPS3?
Our response: We agree with the reviewer on the fact that a bulk crystal of NiPS3 cannot be called a 2D magnet. However, the phenomenon reported here for bulk crystals should be valid for few layer flakes as long as the magnetic order is preserved. This is due to the fact that d-d transitions are localized transitions, and their properties do not depend on the interlayer coupling.
As the magnetic order gets suppressed in the monolayer limit (Ref. 11), the phonon replicas are not expected in the monolayer limit. In addition, phonon replicas have been shown to be existent down to bilayer limit in another publication (see Nature Nanotechnology 16, 655-660, 2021) published during the preparation of this manuscript.

Reviewer #2 (Remarks to the Author):
The manuscript by Ergecen etc., reports the measurements of broadband transient absorption in NiPS3, revealing a phonon replica state. Further coherent phonon spectroscopy shows a phonon oscillation of the same energy as the phonon replica, which only exists within the spectral energy range that corresponds to d-d excitations. The authors conclude that the d-d excitations as the origin of the phonon replica, and deduce electron-phonon coupling strength. The experimental results, like phonon replica and oscillations in coherent phonon spectrum, are clear and convincing. My biggest scientific concern is the lack of temperature dependence in the coherent phonon spectrum, which seems to contradict the behavior of phonon replica and raise more questions in the authors' data interpretation (explained in my detailed comments). Before addressing this main issue, I cannot recommend its publication yet in Nature Communications. Below is my detailed comments: 1. Main concern: Temperature dependence The phonon replica states in the transient absorption spectrum shows the expected temperature dependence, which disappears above TN, because of the symmetry change across the TN as proposed by the authors. But in the coherent phonon spectroscopy (Fig. S6), the oscillations, which is attributed to the same phonon mode (~7THz) as in the phonon replica state, is almost identical as a function of temperature. If they correspond to the same phonon mode, what is the reason for this distinctly different temperature dependence?
The author also mentions that this 7 THz mode agrees with previously assigned Raman modes in Ref 11. In Ref 11, the Ag Raman modes of similar energy all show an obvious temperature dependence, which will agree with the observation in transient absorption spectrum but not the coherent phonon spectrum.
Our response: We thank the referee for examining our work thoroughly and providing insightful comments and inputs.
First, we would like to point out that the Fourier transform of the replica signal is broad and spans an energy interval between 25 meV to 35 meV (Supp. Fig. S5). This range includes both the A1g phonon mode with 253 cm-1 wavenumber (31 meV) and Eg phonon mode with 225 cm-1 wavenumber (28 meV). To pinpoint the phonon mode that is responsible for replica formation, we perform energy-resolved coherent phonon spectroscopy. The most pronounced phonon mode observed in our coherent phonon spectroscopy is the A1g phonon mode with 253 cm-1 wavenumber. In the Raman spectrum, this corresponds to the mode denoted as P5 in Ref. 11. The energy of the A1g phonon mode with 253 cm-1 wavenumber (31 meV) matches the spectral distance between the replicas. On the other hand, the phonon mode at 225 cm-1 (denoted as P4 in Ref. 11), which has negligible spectral weight in Raman measurements, does not appear in our coherent phonon spectroscopy, as shown in Supp. Fig. S7. Therefore, the phonon mode at 225 cm-1 has negligible coupling to the d-d transitions and cannot be responsible for replica formation.
We would like to point out that the temperature dependence of both P4 and P5 phonon amplitudes in Ref. 11, is not correlated with the magnetic ordering temperature. Furthermore, both of them exist above and below the magnetic ordering temperature. This observation indicates that the Raman amplitude for 253 cm-1 phonon mode is not directly influenced by the magnetic order. As we have mentioned in the main text, the coupling between the d-d transitions and 253 cm-1 phonon mode exists above the magnetic ordering temperature, but they are not optically active because of local inversion symmetry. At low temperatures, d-d transitions become optically active because of local inversion symmetry breaking arising from the magnetic order.
As the referee alluded to, the Raman amplitudes of both P4 (225 cm-1) and P5 (253 cm-1) phonon modes show temperature dependence, and this dependence is not observed in our coherent phonon spectroscopy data. Even though both Raman and coherent phonon spectroscopy are sensitive to phonon modes, the excitation and detection of the phonon modes are completely different for these spectroscopy modalities. In Raman, a single photon inelastically scatters and spontaneously creates a single phonon excitation through a virtual transition, highly detuned from equilibrium excited states. Thus, the Raman amplitude depends on the structure of the equilibrium excited state (detuning etc.) and the thermal occupation of phonon modes, which can change as a function of temperature.
On the other hand, in coherent phonon spectroscopy, a probe pulse samples the real time phase coherent phonon oscillations following a pump excitation. Unlike Raman scattering, these coherent phonon oscillations are launched by impulsive Raman scattering and displacive excitation of coherent phonons (DECP). In the case of DECP, the amplitude of the phonon mode is proportional to the phonon displacement induced by the pump pulse, which is determined by the difference between equilibrium and nonequilibrium lattice positions. Our coherent phonon spectroscopy data for NiPS3 implies that the excitation amplitude of the P5 phonon mode does not change with temperature. We think that the discrepancy in temperature dependences of the Raman and coherent phonon spectroscopy is not unexpected, as they can launch and probe phonon modes differently.
To further clarify the differences between Raman and coherent phonon spectroscopy, we would like to compare and contrast two studies (one Raman and one coherent phonon spectroscopy) performed on isostructural & isoelectronic van der Waals magnets CrSiTe3 & CrGeTe3. For these compounds, the A1g phonon mode at 136 cm-1 does not show any temperature dependence in Raman (Fig. 4 -https://arxiv.org/pdf/1604.08745.pdf), whereas coherent phonon spectroscopy (https://arxiv.org/pdf/1910.06376.pdf) shows a dramatic change in amplitude because of a spin dependent phonon excitation mechanism. We hope that this example clears up the differences between coherent phonon spectroscopy and Raman measurements.
We have added a supplementary note (Supp. Note 10) that describes the differences between Raman and coherent phonon spectroscopy more clearly, and made the following changes in the main text: 1) [They manifest themselves… -Line 38] -We clarified the frequency of the phonon mode that participates in the replica formation. Minor comments/concern: 2. Page 3, 2nd line '… with strong vibronic coupling [10]'. I'm not sure if this ref 10 is the best for this purpose. Our response: This was a mistake. We have corrected it by replacing the reference with the right one. We thank the referee for pointing this out. 3. Compared to Fig. 2, it's hard to resolve the phono replicas in Fig.S4. Is it just because of the color scaling difference? Our response: We agree with the referee. In Fig. S4, it is harder to resolve the phonon replicas. The data in Fig. 2, was taken by fixing the delay time (at 2 ps), and averaging for a longer time. On the other hand, the data in Fig. S4 had been averaged less for the sake of time, since we are examining the time dependence as well. 4. When extracting the phonon replica energy and oscillation frequency, please also add the error bar the number. It seems to have quite some energy distribution. Our response: We thank the referee for bringing this to our attention. Following this comment, we made the following changes: 1) [Supp. Note 9] -we have added the error bars to the extracted values of replica energy and oscillation frequency for two different undressed d-d transition lineshapes, Gaussian and Lorentzian. 2) [The energy spacing between -Line 66] -We rephrase this sentence to better describe the energy distribution of the replicas seen in our transient absorption measurements.
5. In Fig. 2b, the arrows indicating energy spacing can be misleading. It should have corresponded to the length start/end at the center of the replica. The size of the arrow (compared to x scale) is also smaller than the extracted 28.5meV. Our response: We have corrected this in our revised version of the manuscript. The arrows were for demonstration purposes only. We removed them in our revised manuscript. Thanks for pointing this out.

Reviewer #3 (Remarks to the Author):
Report for "Magnetically brightened dark electron-phonon bound states in a van der Waals antiferromagnet".
In this manuscript, Ergecen and Ilyas et al revealed a unique electron-phonon bound state below the Néel transition in NiPS3, and more importantly, identified the electronic origin as localized d-d excitations and the phonon source as A1g phonon at ~ 7THz, using transient absorption spectroscopy with a dynamic range. Van der Waals (vdW) magnets emerge as a flourish platform that hosts strong coupling between multiple degrees of freedom and realizes novel bound states down to the two-dimension (2D) limit. Yet, it is highly challenging to detect the presence of these bound states and more so to resolve the composition of them. This manuscript has shown, for the first time in vdW magnets, the unambiguous identification of electronic and phononic origins of an electron-phonon bound state with the strongest electron-phonon coupling known so far. This manuscript is of high interest and high quality. I would recommend it for publication at Nature Communications, after addressing the following questions.
1. Could the authors please provide an equilibrium absorption spectrum? How does this equilibrium spectrum compare with the transient absorption spectra? Are the fine oscillations for the electron phonon coupled states resolvable in the equilibrium spectrum? Our response: We are thankful to the referee for carefully reading our work, and appreciating the importance of our findings.
The equilibrium absorption spectrum of NiPS3 has been reported in a recent publication (Phys. Rev. Lett. 120, 136402 (2018) -Fig. S6). The equilibrium absorption spectrum cannot resolve the fine spectral oscillations. This fact highlights the power of transient absorption measurements, which have a high dynamic range and allow us to observe faint spectral details.
In addition, equilibrium linear dichroism (birefringence) spectroscopy measurements have also observed the phonon replicas (see Nature Nanotechnology 16, 655-660, 2021) published during the preparation of this manuscript. Because NiPS3 exhibits magnetically induced birefringence, equilibrium linear dichroism (birefringence) spectroscopy is sensitive to the magnetic order induced spectral features with high dynamic range.
2. What is the role of demagnetization with the high intensity pump in detecting the phonon replicas in the transient absorption spectroscopy? Why is the transient reflectivity probed 2 picoseconds after the pump? What is the temporal dynamics of the phonon replicas?
Our response: The high intensity pump pulse heats up the electronic system by efficiently generating electron-hole pairs. This in turn melts the magnetic order and "washes out" any spectral features pertinent to magnetic order. We measure the difference between this nonequilibrium absorption spectrum (at 2 ps time delay) and equilibrium one (at negative time delay). This is what allows us to measure small spectral signals, which cannot be detected with conventional equilibrium absorption methods.
The 2 ps time delay is not a special point. We actually have examined the temporal dynamics of phonon replicas up to 100 ps delay times (see Fig. S4) and have observed that phonon replicas survive up to this delay times with no visible change in oscillation amplitudes. This indicates that the spectral oscillations are indeed equilibrium phenomena.
3. When fitting the spectrum with a sum of Gaussians weighted by a Poisson distribution to extract the Huang-Rhys factor g, there are a few questions that I am interested in: a. Why do the authors choose Gaussian, instead of Lorentzian? Our response: The model we use to obtain the Huang-Rhys factor g takes the undressed d-d transition energy as a free parameter. In our fitting procedures, we both used Gaussian and Lorentzian distributions for the undressed d-d transition lineshapes. The selection of the lineshape function does not significantly change the extracted Huang-Rhys factor and the phonon frequency. Therefore, none of the conclusions are affected by the selection of the lineshape function.
We have added a subsection into Supp. Note 9 that shows the fitting results obtained using a Lorentzian lineshape.
b. How do the authors determine where the zeroth order of the Poisson distribution is? Our response: In our fits, the model we use takes the bare d-d transition energy as a free parameter. The fit outputs the zeroth order of the Poisson distribution.
c. Is there a broad background in the absorption spectrum? What happens if fitting with a broad background plus a sum of Gaussians weighted by a Poisson distribution? Our response: The model used in this paper fits the d-d transition region very well without any need of a background.
d. The linewidth of individual Gaussians should be the convolution between the phonon linewidth and the d-d excitation linewidth. However, the d-d excitation linewidth, reported in PRL 120, 136402 (2018), seems to be much larger than the Gaussian linewidth here. Any reasons?