Time and momentum-resolved tunneling spectroscopy of pump-driven non-thermal excitations in Mott insulators

TL;DR

This paper introduces a computational tunneling spectroscopy method to analyze non-thermal excitations in Mott insulators, revealing detailed momentum-resolved spectral features post-quench.

Contribution

The authors develop a novel extended probe technique for momentum-resolved spectral analysis in correlated systems using a tunneling approach.

Findings

Identification of in-gap sub-bands associated with non-thermal states

Detection of excitons and anti-bound states at high energies

Observation of negligible relaxation within the studied timescales

Abstract

We present a computational technique to calculate time and momentum resolved non-equilibrium spectral density of correlated systems using a tunneling approach akin scanning tunneling spectroscopy. The important difference is that our probe is extended, basically a copy of the sample, allowing one to extract the momentum information of the excitations. We illustrate the method by measuring the spectrum of a Mott-insulating extended Hubbard chain after a sudden quench with the aid of time-dependent density matrix renormalization group (tDMRG) calculations. We demonstrate that the system realizes a non-thermal state that is an admixture of spin and charge density wave states, with corresponding signatures that are recognizable as in-gap sub-bands. In particular, we identify a band of excitons and one of stable anti-bound states at high energies that gains enhanced visibility after the pump. We do not appreciate noticeable relaxation within the time-scales considered, which is attributed to the lack of decay channels due to spin-charge separation. These ideas can be readily applied to study transient dynamics and spectral signatures of correlation-driven non-equilibrium processes.

Equations

$$H_0=-J\sum _{i=1,\sigma }^{L-1}\left(c_i^{†}{c}_{i+1}+h.c.\right)=\sum _k{\omega }_k{c}_k^{†}{c}_k,$$

$$H_{\mathrm{probe}}=V_g\sum _{i=1}^L{d}_i^{†}{d}_i=V_g\sum _k{d}_k^{†}{d}_k,$$

$$H_{\mathrm{tunnel}}=J^{\prime }\sum _{i=1}^L\left(d_i^{†}{c}_i+h.c.\right)=J^{\prime }\sum _k\left(d_k^{†}{c}_k+h.c.\right).$$

$$H=\left[\begin{array}{cc}{\omega }_k& -J^{\prime }\\ -J^{\prime }& V_g\end{array}\right]$$

$$E_{\pm }(k)=\frac{V_g+{\omega }_k}{2}\pm \sqrt{{\left(\frac{V_g-{\omega }_k}{2}\right)}^2+J^2},$$

$$⟨n_d(k,t)⟩={\int }_0^{t}d{t}_1{\int }_0^{t}d{t}_2⟨V_k(t_1)n_d(k)V_k(t_2)⟩,$$

$$⟨n_d⟩$$

$$H$$

$$n_{d\sigma }(k)=\frac{1}{L}\sum _{j,l=1}^L{e}^{ik(j-l)}⟨d_{j\sigma }^{†}{d}_{l\sigma }⟩.$$

Full text

I Introduction

With the advent of new powerful light sources, experimentalists can shake the excitations of a system and probe states present in the spectrum that are not accessible via finite-temperature measurements. By means of ultrafast light pulses, electrons can be excited above any intrinsic energy scale, and the competition between different degrees of freedom can be manipulated. Cavalieri et al. (2007); Orenstein (2012); Chollet et al. (2005); Corkum and Krausz (2007); Onda et al. (2008). The resulting non-thermal states after photoexcitation often contain coexisting orders that are not usually present in the ground or thermal statesOrenstein (2012); Kampfrath et al. (2013). This new knob can be used to stabilize “hidden” phases that reside at higher energies, such as superconductivityFausti et al. (2011), and to induce or disrupt charge, magnetic, or orbital order Polli et al. (2007); Ehrke et al. (2011); Zhang et al. (2016); Casals et al. (2016); Chollet et al. (2005); Onda et al. (2008); Okamoto et al. (2007).

Time-resolved femtosecond photoemission spectroscopy has been one of the most a used techniques to monitor in real time and with atomic resolution the ultrafast quasiparticle dynamics in correlated-electron materials Plummer (1997); Ramakrishna et al. (2001); Domcke (1991); Sebastian and Tachiya (2006). The experimental protocol Cavalieri et al. (2007); Smallwood et al. (2016) starts with an intense pulse of radiation that ‘pumps’ the system into a highly excited non-equilibrium state. After a variable time delay, the system is subject to a weak probe pulse of higher energy photons, ejecting photoelectrons which are detected with energy (and angle) resolution. By means of this powerful tool, one can peek into the different decay mechanisms taking place, and experimentally unveil the complex and rich interplay between charge, spin, orbital and vibrational degrees of freedom.

Notwithstanding, theoretically reproducing time- and angle-resolved photoemission spectra is computationally challenging and expensive. It can be numerically carried out only in small systems, as it requires the full knowledge of the eigenstates and the calculation of a two-time correlator Freericks et al. (2009); Shao et al. (2016). In the equilibrium steady state, approximations can be made by using the single-particle Green’s function, but all information about transient and the actual decay mechanisms during the relaxation process is lost. In a non-thermal state far from equilibrium, the imaginary part of the equilibrium retarded Green’s function $G(\omega)$ is not guaranteed to be positive and does not yield meaningful information about the orbital occupation (it is not a density of states).

We hereby propose a different approach to investigate these quantities using a tunneling technique. We focus on a geometry that was first suggested in Ref.Carpentier et al., 2002, and later realized experimentally in Refs.Auslaender et al., 2002, 2005 for conducting momentum-resolved tunneling spectroscopies on one-dimensional (1D) systems. Unlike scanning tunneling spectroscopy, where the probe yields only local informationCohen et al. (2014), an extended one-dimensional wire can provide momentum resolution. Electrons can tunnel from the sample into the one-dimensional non-interacting lead that is placed parallel to it. Since this occurs in the transverse direction, momentum conservation along the probe direction is ensured. A gate voltage $V_{g}$ is applied to the probe wire and energy conservation implies that only electrons with energy $\omega=V_{g}$ can tunnel. Momentum resolution is achieved through the application of a magnetic field perpendicular to the plane of the sample and probe wires. A similar scheme was recently proposed for performing momentum-resolved spectroscopies on cold atomic systems: instead of a voltage, an RF field or the shaking of the lattice can yield transitions at a target frequencyKantian et al. (2015); Bohrdt et al. (2018). In this variation, as particles tunnel to the second channel, momentum is mapped via time of flight.

In section II we describe in detail the implementation and illustrate with simple examples. In section III we present results for an interacting system –the extended Hubbard model– after a quench, and we close with a summary and discussion.

II Method

We propose to computationally carry out a hybrid method combining ideas from the aforementioned setups: after the system has been photoexcited, we allow for electrons to tunnel into an empty parallel wire which has been set at a given gate voltage, as shown in Fig.1(a). Only electrons at a particular energy $V_{g}$ can tunnel, and we can then access the occupation of each state with momentum resolution by simply calculating the momentum distribution function of the probe wire.

II.1 Non-interacting fermions

We illustrate this idea with the simple example of non-interacting fermions, whose Hamiltonian reads

[TABLE]

where $c^{\dagger}_{i}$ and $c_{i}$ are the usual creation and annihilation fermion operators (we ignore the spin index for now) and $w_{k}=-2J\cos{k}$. We take the inter-atomic distance as unity and we express all energies in units of the hopping parameter $J$ (the symbol “$t$” will be reserved to represent time, which will be expressed in units of $1/J$).

A second “probe” chain is included as

[TABLE]

where we distinguish the operators $d^{\dagger}$ and $d$ acting on the probe. Notice that there is no hopping nor interactions along the probe chain: it consists of isolated empty orbitals with a gate voltage (or chemical potential) $V_{g}$. At time $t=0$ the system is in the ground state of the physical chain at a fixed given density, while the probe chain is empty (this is ensured by initially setting $V_{g}$ to a very large positive value). Then, both chains are connected by means of a tunneling term:

[TABLE]

Putting together Eqs.(1),(2) and (3), the full problem becomes the sum of $L$ independent tunneling terms, which can be readily solved. For simplicity, we look at an eigenstate at temperature $T=0$, in which a single particle orbital with momentum $k$ is either empty or occupied. The single particle states are $|n(k),n_{d}(k)\rangle={|1,0\rangle;|0,1\rangle}$ where $1$ and [math] represent the occupancy of the physical orbital or the probe orbital with momentum $k$. The Hamiltonian for the two level system is :

[TABLE]

with ground state energies

[TABLE]

Starting from a initial state $|1,0\rangle$, the probability that a fermion is transferred to the corresponding empty probe state $k$ at time $t$ is simply $n_{d}(t)=|\langle t|0,1\rangle|^{2}=4A^{2}\sin^{2}{\left(\frac{E_{+}-E_{-}}{2}t\right)}$, with $A=J^{\prime}W/(J^{\prime}+W)^{2},W=\omega_{k}-E_{+}$. This function oscillates in time with a period $\tau=\pi/J^{\prime}$ for $V_{g}=\omega_{k}$. In order to maximize the “visibility” one needs to measure the density of the probe state $k$ at time $t_{\mathrm{max}}=\pi/2J^{\prime}$. As a function of $V_{g}$, the probability is peaked at $V_{g}=\omega_{k}$ (with smaller satellite peaks), and its width gets narrower as $J^{\prime}\rightarrow 0$, or $t_{\mathrm{max}}\rightarrow\infty$, as shown in Fig.2. This is nothing else but Fermi’s golden rule and a manifestation of the uncertainty principle: to obtain sharper resolution in energy, one needs to choose a small coupling $J^{\prime}$ and measure at very long times.

II.2 General formulation away from equilibrium

We now consider the full many-body Hamiltonian $H$ and a generic initial state $\rho_{0}=\rho\otimes|0\rangle\langle 0|$. The density matrix $\rho$ represents the state of the many-body system: $\rho=\sum_{n,n^{\prime}}a_{n^{\prime}}^{*}a_{n}|n^{\prime}\rangle\langle n|$, where the states $|n\rangle$ are the eigenstates of $H$. We point out that this is a general scenario, in which the system may have been driven away from equilibrium by an external perturbation or a quench and $H$ is the final Hamiltonian; the equilibrium case is simply recovered by taking a diagonal density matrix. We assume that the measurement process starts suddenly some time after the perturbation, which for simplicity of notation we label as $t=0$, and the system evolves thereafter under the action of a time-independent Hamiltonian (i.e.,any time-dependence in the Hamiltonian is “frozen” at time $t=0$). Following closely the discussion in Ref.Kantian et al., 2015 we find that in second order of perturbation, the occupation of the state $k$ of the probe system is given as (we ignore the spin index for now):

[TABLE]

where $V_{k}=-J^{\prime}\left(d^{\dagger}_{k}c_{k}+\mathrm{h.c.}\right)$. The averages are with respect to the initial state $\langle\cdots\rangle=\mathrm{Tr}(\cdots\rho_{0})$ and we work in the interaction picture $O(t)=e^{i(H+V_{g}n_{d})t}Oe^{-i(H+V_{g}n_{d}t)t}$ (from now on we ignore the subindex $k$ for convenience). Since the probe orbital is initially empty, the only term surviving in this expression is:

[TABLE]

Moreover, noticing that the initial state is a product state, we can readily evaluate the contribution of the probe orbital to this expression: $\mathrm{Tr}_{probe}(d(t_{1})n_{d}d^{\dagger}(t_{2})|0\rangle\langle 0|)=e^{iV_{g}(t_{2}-t_{1})}$. Explicitly, Eq.(5) becomes:

[TABLE]

with $\omega_{nm}=E_{n}(N)-E_{m}(N-1)$, $N$ being the number of particles. Hence, the resulting occupation will be peaked at the gate voltages corresponding to the allowable transitions, $\langle m(N-1)|c_{k}|n(N)\rangle$ , weighed by the initial occupation of the eigenstates. We notice that the argument in the integral is just the lesser Green’s function $G^{<}(t_{1},t_{2})=\langle c^{\dagger}(t_{1})c(t_{2})\rangle$ and this equation is identical to the one derived in Ref.Freericks et al., 2009 to describe a time-resolved photoemission experiment.

For large $t$, the quantity $n_{d}(k,t)$ converges to a sum of Dirac deltas and yields an expression proportional to the system’s spectral function. Clearly, at long times the electron will be reflected and tunnel back to the system so, in reality, to improve the energy resolution one needs to pick $J^{\prime}$ small. In general, we take as a rule of thumb $t_{\mathrm{max}}=\pi/2J^{\prime}$ in all cases. We point out that a time-dependent tunneling term could also be considered, which translates into the introduction of an envelope function, as done in Ref.Freericks et al., 2009.

III Results

We now demonstrate an application of this scheme to explore competing orders and excitations in one-dimensional correlated materials. It is known that in 1D systems, the band edge singularity could give rise to a high-differential optical gain, with potential applications such as light-emitting diodes, lasers, sensors, and molecular switches Brown et al. (1995); Burroughes et al. (1990); Dodabalapur et al. (1995a, b); Hide et al. (1996); Yang and Swager (1998); Schmitz et al. (2001); Nitzan and Ratner (2003). There is great deal of interest in the optical properties of 1D materials in the presence of correlations, when a gap arises as a result of electronic interactions. Moreover, the emergence of excitonic excitations, has been subject of attention of a number of theoreticalJeckelmann et al. (2000); Tsutsui et al. (2000); Essler et al. (2001); Jeckelmann (2003); Gallagher and Mazumdar (1997); Barford (2002); Gebhard et al. (1997); Kancharla and Bolech (2001); Mizuno et al. (2000); Glocke et al. (2007); Matsueda et al. (2004); Lu et al. (2015); Al-Hassanieh et al. (2008); Dias da Silva et al. (2010); Rincón et al. (2014)and experimental Ono et al. (2005); Schlappa et al. (2012) works.

The minimal model to study correlated polymers is the so-called “$U-V$” extended Hubbard model:

[TABLE]

Here, $c^{\dagger}_{i\sigma}$ creates an electron of spin $\sigma$ on the $i^{\rm th}$ site along a chain of length $L$. The on-site and nearest-neighbor Coulomb repulsion are parametrized by $U$ and $V$, respectively.

The physics of one-dimensional strongly correlated fermionic systems can generally be described in terms of Luttinger liquid theory. In a Luttinger liquid (LL) Haldane (1981); Gogolin et al. (1998); Giamarchi (2004), the natural excitations are collective density fluctuations, that carry either spin (“spinons”), or charge (“holons”). This leads to the spin-charge separation picture, in which a fermion injected into the system breaks down into excitations, each with a characteristic energy scale and velocity (one for the charge, one for the spin). Spin-charge separation acts as a constraint for the dynamics of the system, that cannot relax to a thermal state after a quench or non-equilibrium situation. The lack of thermalization implies that it might be possible to ‘trap’ the system in an excited state for very long times.

As a proof of concept we conduct a numerical experiment using the time-dependent density matrix renormalization group method (tDMRG) White and Feiguin (2004); Daley et al. (2004); Feiguin (2011); Paeckel et al. (2019) on chains of length $L=32$ and with parameters $U=20$ and $V=5$. This choice may seem exaggerated, but is justified: it will provide us with a large Mott gap $\Delta\sim U-V$, and allow us to resolve any features that may appear inside the gap with more detail and well separated from the bands. The ground state of the system at half-filling is a Mott insulator with dominant power-law decaying quasi-long-range antiferromagnetic order, or SDW phase Nakamura (2000); Jeckelmann (2002); Sandvik et al. (2003); Tsuchiizu and Furusaki (2004). The optical conductivity and Raman spectrum reveal the existence of sharp excitonic peaks with a weak continuous band of free excitations of width $\sim 8J$ Jeckelmann et al. (2000); Tsutsui et al. (2000); Essler et al. (2001); Jeckelmann (2003); Gallagher and Mazumdar (1997); Gebhard et al. (1997); Kancharla and Bolech (2001). However, these optical excitations are not present in the spectrum, shown in Fig.3(a) as a reference, also obtained using tDMRG with $m=600$ states. The lower and upper Hubbard bands are well separated from each other by a wide Mott gap, and no remarkable features are observed, besides the characteristic holon and spinon dispersions.

To avoid considerations concerning pulse shape, frequency, and length, we simplify the discussion to the case of a quench, in which the system is prepared in the Mott insulating ground state of a system of $N=L$ electrons with $U_{0}=2,V_{0}=0$, and the interactions are suddenly changed to $U=20,V=5$. As a consequence, the final state will be a superposition of eigenstates that will exhibit a large number free holes and doublons, as well as excitons, occupying broad range of energies. In Fig.4 we show results obtained using tunneling spectroscopy right after the quench. The probe is connected to the chain at time $t_{\mathrm{wait}}$ after the quench, and we plot the momentum distribution function of the probe chain at time $t_{\mathrm{wait}}+t_{\mathrm{probe}}$ as a function of momentum and gate voltage:

[TABLE]

Notice that we use open boundary conditions throughout, which translates into some uncertainty in momentum. We scanned $V_{g}$ in steps of 0.2, implying 175 independent tDMRG simulations for each value of $t_{\mathrm{wait}}$. We took $J^{\prime}=0.2$ and used $m=200$ DMRG states, which yields a truncation error of the order of $10^{-4}$ in the worse cases.

As shown in Fig.4, besides some sharper and better defined features, we are not able to resolve a noticeable difference between the measurements right after the quench and at $t_{\mathrm{wait}}=5$. This is also reflected in the integrated weight over momenta, displayed in Fig.5: panel (a) illustrates how the visibility improves as a function of $t_{\mathrm{probe}}$ (see animations in the supplementary material), while in (b) we compare the two waiting times. In this case, we are able to resolve some minor differences that stem from the relaxation of excitations within the lower Hubbard band, indicating the lack of available channels for non-radiative decay or recombination. It is possible that these excitations cannot decay due to the energy mismatch between the bandwidth $W$ and the interactions, or thermalization occurs in timescales that far exceed the simulation time. This bottleneck exists already in higher dimensions Sensarma et al. (2010); Eckstein and Werner (2011); Lenarčič and Prelovšek (2013); Eckstein and Werner (2016). If the bandwidth is small, the number of available decay channels gets suppressed. However, in our case the in-gap states are not too far from each other, nor from the lower Hubbard band. In higher dimensions it was observed that the spin excitations are highly relevant for thermalization. It is possible that spin-charge separation, which is more dramatic at large values of $U$, and the flat spinon dispersion for large values of $U$ do not allow for a wide range of energy and momenta for scattering.

The spectrum is very well resolved and displays many non-trivial features that are not present neither in the zero temperature spectrum nor the optical conductivity. In order to account for these results, we first assume the possibility that the system is in a thermal state. We have calculated the spectra for a wide range of temperature scales and have found that the final state after the quench does not correspond to a thermal distribution. For illustration, we display finite-temperature tDMRGFeiguin and White (2005) results at $T=2.5J$ in Fig.3(b). The first remarkable and most obvious feature of the spectrum is recognizable in the lower Hubbard band, which displays a dispersion rather resembling a tight-binding band of spinless fermions than the usual characteristics of fractionalized excitations seen in panel (a). This is actually expected, since in this regime the spin is completely incoherent (We refer the reader to Refs.Cheianov and Zvonarev, 2004; Cheianov et al., 2005; Abendschein and Assaad, 2006; Fiete, 2007; Halperin, 2007; Feiguin and Fiete, 2010, 2011; Soltanieh-ha and Feiguin, 2014; Nocera et al., 2018 for a discussion of the finite-temperature spectra of 1D correlated systems). Moreover, we distinguish a distribution of spectral weight inside the gap due to the correlated nature of the problem Nocera et al. (2018), a phenomenon that has been experimentally observed in the photoemission spectrum of the single chain Mott insulators Sr2CuO2Kidd et al. (2008) and Na0.96V2O5 Kobayashi et al. (1999). On the other hand, the tunneling spectrum displays a quite large spectral weight inside the gap and in the upper Hubbard band, implying that if we had to assign a temperature to the system after the quench, it would have to be larger than the Mott gap. However, unlike the finite temperature case, the spinon and holon bands remain coherent.

In order to make sense of the unexpected features in the tunneling results, we carry out a similar simulation using exact diagonalization on a chain with $L=4$ sites with a parallel chain as a probe. The complexity of the problem is similar to that of a $2\times 4$ Hubbard ladder with 4 electrons. Even though it is a small system and is likely very affected by boundary effects, it provides valuable intuition to interpret the tDMRG results. Following a similar protocol, we first resolve the tunneling spectrum, shown in Fig.6(a). Since we have access to all eigenstates and eigenvalues, we calculate all possible single particle excitation energies as $\omega_{nm}=E_{n}(N=L-1,S^{z}=1)-E_{m}(N=L,S^{z}=0)$, some of which are shown in the plot with different colors. The final state is predominantly a superposition of the ground state –which has dominant SDW correlations– and two excited states, labelled $|m=3\rangle$ and $|m=5\rangle$ in Fig.6(b), that display CDW correlations, as shown in panel Fig.6(c). This enhancement of the charge order was previously observed in Ref.Lu et al., 2012 under the action of a driving field. We focus on the dominant features of the spectrum, namely, the flat bands at energy $\omega\sim-5$ and $\omega\sim-15$, and the in-gap spectral weight at energies between $\omega=5$ and $\omega=7$. The first one corresponds to breaking a holon-doublon pair on top of $|m=3\rangle$, while the in-gap weight corresponds to excitations on top of $|m=5\rangle$. The flat band at high energies below the Fermi level is an excitation on top of the ground state that acquires an enhanced spectral weight after the pump. This high energy feature has been overlooked in prior studies of the model due to its very weak spectral signatures at zero-temperature, and indicates the presence of stable anti-bound states outside of the continuum.

IV Conclusions

To summarize, we have introduced a computational tunneling approach that allows one to access the time and momentum resolved spectrum of strongly correlated systems away from equilibrium, which previously could only be obtained from small systems with exact diagonalization. The formulation is general and does not depend on how the system is driven out of equilibrium. We have applied the method to study the dynamics of Mott insulating Hubbard chains after a quench and have been able to identify features in the spectrum corresponding to an admixture of SDW and CDW states, with a band of doublon-holon excitons and high-energy anti-bound states. This extremely powerful technique can be readily extended to arbitrary models under a variety of scenarios, giving access to transient dynamics and the ability to identify correlation-driven non-equilibrium processes behind pump-driven phase transitions and exciton decay and recombination.

Acknowledgements.

We thank A. Nocera for valuable comments. We acknowledge generous computational resources provided by Northeastern University’s Discovery Cluster at the Massachusetts Green High Performance Computing Center (MGHPCC). KZ is supported by a Faculty of the Future fellowship of the Schlumberger Foundation. AEF acknowledges the U.S. Department of Energy, Office of Basic Energy Sciences for support under grant No. DE-SC0014407.