Distinct charge and spin recovery dynamics in a photo-excited Mott insulator

With the discovery of high temperature superconductivity in the doped cuprates their parent Mott insulating state, for example La₂CuO₄ [1], has been extensively studied [2]. The Mott state arises due to strong local repulsion in atomic orbitals, which prevents simultaneous occupancy of both spin states, leading to an insulator at half-filling [3, 4, 5]. The localised electrons have an inter-site exchange interaction that promotes antiferromagnetic order in non frustrated lattices.

When the Mott insulator is excited by a laser pulse with appropriately chosen frequency the added energy has an impact on both the charge dynamics and the magnetic order [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. The pump can excite electrons from the lower Hubbard band (LHB) to the upper Hubbard band (UHB), which in real space means creation of double occupancy [22, 23, 24]. The effect is twofold: (i) the electrons excited to the UHB, and the ‘holes’ in the LHB, act as mobile carriers, leading to metallic response in the optical conductivity, and (ii) the double occupancy suppresses the magnetic moments, destroys their spatial correlation, and leads to suppression of magnetic long range order. In effect, despite being at half filling, one obtains a transient metallic state coexisting with (small) local moments, which evolves back towards its reference antiferromagnetic Mott state. It is in this context that Sr₂IrO₄ has thrown up several puzzles about charge and spin dynamics out of equilibrium.

A photo-excited state was realised in Sr₂IrO₄, a structural analog of La₂CuO₄. At equilibrium Sr₂IrO₄ is a spin-orbit coupled layered antiferromagnetic Mott insulator (AFMI) with an interaction split $J_{eff}=1/2$ band [25, 26, 27, 28]. The AF $T_{c}$ is $\sim 240$K. Laser pumping leads to dramatic changes in its magnetic and electronic properties. The two most significant observations in our reading are: (i) The rapid loss of 3D magnetic order and gain in reflectivity in response to a pump - and then a quick suppression of the reflectivity, within 1-2 picoseconds (ps) but very slow recovery of order (touching 1000 ps) [29]. This is the ‘two timescale’ issue. (ii) The persistence of low energy spectral weight in the optical conductivity even after a seeming steady state is quickly reached, suggesting a population of excited high energy electrons - the holon-doublon (HD) plasma - at long times [30]. The time dependence of planar magnetic fluctuations [29], oddly, reflects both short and long timescales! In the current interpretation [29] the in plane magnetic order and electron physics in the Mott insulator both recover quickly and the 3D recovery is delayed due to weak interlayer coupling.

We suggest that the fascinating data revealed by the experiments has another - very different - explanation. This has remained out of reach because theories of nonequilibrium phenomena need to consider the excited electronic population and also the spatial dynamics, e.g, domain growth effects, on large spatial scales. Current tools, e.g, exact diagonalisation [31, 32], dynamical mean field theory (DMFT) [33, 34, 35, 36, 37, 38]. and density matrix renormalisation group [39, 40], are either size limited or unable to access the timescale needed.

To capture the electron correlation non perturbatively and access spatio-temporal dynamics we write an equation of motion for the one body density operator $\rho_{ij}^{\sigma\sigma^{\prime}}=c^{\dagger}_{i\sigma}c_{j\sigma^{\prime}}$ in the half-filled square lattice Hubbard model and close the hierarchy by factorising the resulting four operator terms in the magnetic channel. This is ‘mean field dynamics’ (MFD) although it is in terms of the matrix $\rho_{ij}^{\sigma\sigma^{\prime}}$ rather than a local object like magnetisation. This allows access to size $\sim 20\times 20$ and time upto $10^{3}t^{-1}_{hop}$, where $t_{hop}$ is the in-plane hopping. We also construct, benchmark, and use a Langevin dynamics (LD) scheme that allows access to $\sim 60\times 60$ lattices and time upto $10^{4}t^{-1}_{hop}$. We set $t_{hop}=260$meV ($t_{hop}^{-1}\equiv 16$ femtoseconds) and $U=3t_{hop}$ as noted [25] for Sr₂IrO₄. Our results, based on a combination of MFD and LD, are the following.

(i) The pump induced suppression of order and enhancement of optical weight occurs over $\sim 0.3$ps. Post pump, the optical weight reduces to reach its long time asymptote over $\tau_{opt}\sim 2-3$ ps, while magnetic order recovery time $\tau_{ord}$ ranges from a few ps to $100$ps depending on fluence.

(ii) The ‘time resolved’ magnetic fluctuation spectrum $S({\bf q},\omega,t)$ has a recovery time $\tau_{fluc}({\bf q})$ that varies widely, $\tau_{fluc}({\bf q})\sim\tau_{opt}$ when ${\bf q}$ is far from the ordering vector ${\bf Q}=(\pi,\pi)$, and $\tau_{fluc}({\bf q})\sim\tau_{ord}$ when ${\bf q}\rightarrow(\pi,\pi)$.

(iii) While ‘charge recovery’ to a steady state value, seen in $\tau_{opt}$, is quick, a significant UHB electron population survives to long time, sustaining low frequency optical weight. $\tau_{opt}$ and the weight correlates with measurements [30].

(iv) Studying a layered $O(3)$ symmetric model with $J_{\perp}/J_{\parallel}\sim 10^{-3}$ we found that the 3D ‘recovery time’ differs from the 2D value by only a factor $\sim 2$. Key dynamical features of the layered material are two dimensional.

Model: Although Sr₂IrO₄ is a three dimensional multiband system with spin-orbit coupling the essential physics is in the IrO planes, controlled by a single spin-orbit coupled orbital per site subject to a Hubbard interaction. The single band model takes the conventional form: $H=\sum_{\langle ij\rangle,\sigma}t_{ij}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}$. We set the lattice spacing to 1.

We introduce the pump via a pulse of amplitude $E_{0}$, frequency $\Omega_{p}\sim 120$ THz and pulse width $\tau_{p}=100$fs by Peierls coupling. $E_{0}^{2}$ would be proportional to the fluence. The evolution of $\langle\rho_{ij}^{\sigma\sigma^{\prime}}\rangle$ involves solution of the $4N^{2}$ MFD dynamical equations [41], where $N$ is the number of sites. The form of this equation, its derivation, and its implementation is discussed in Supplement A [42], and primary indicators shown in Supplement B. The local magnetisation can be calculated as ${\vec{m}}_{i}=\sum_{\sigma,\sigma^{\prime}}\vec{\tau}_{\sigma\sigma^{\prime}}\langle\rho_{ii}^{\sigma\sigma^{\prime}}\rangle$, where $\tau$ are the Pauli matrices.

While MFD reveals key features of ‘suppression-recovery’ dynamics in response to the pump, it has limitations in terms of accessible size and time that prevent study of the timescale separation observed in experiments. Electronic properties, like the upper Hubbard band occupancy, stabilize to a pump-dependent constant value on a short timescale $t\sim\tau_{opt}\sim 2$ ps, as confirmed within MFD to $t\sim 10\tau_{opt}$. Beyond a few ps, the population of electrons excited to the upper Hubbard band stabilizes to a finite time-independent value, akin to having a steady finite ’electronic temperature’ $T_{el}$ at long times. This constancy is then used in Langevin calculations on large lattices for times up to $(50-60)\tau_{opt}$. The detailed characterization of $T_{el}(t,E_{0})$ is provided in Supplement B, where $T_{el}$ serves as a crucial input for modeling magnetic dynamics below.

To address the magnetic dynamics with high spatial resolution and long times we write a Langevin equation directly for the ${\vec{m}}_{i}$. In such an equation the ${\vec{m}}_{i}$ are subject to a torque arising from the electrons, a damping, and a ‘thermal’ noise at the bath temperature $T_{b}$. The pump excited electrons follow a Fermi distribution with temperature $T_{el}$ inferred from MFD. The equation takes the form:

$F_{SF}$ is the free energy of electrons at a temperature $T_{el}$ in the spin background $\{{\vec{m}}_{i}\}$, $\gamma$ is a dissipation constant extracted from MFD, and $\eta_{i}$ is thermal noise with $\eta$ and $\gamma$ satisfying the fluctuation-dissipation theorem at temperature $T_{b}$. The free energy $F_{SF}$ arises from the Hubbard model, $H_{SF}$, written in the spin-fermion language:

${\vec{s}}_{i}=c^{\dagger}_{i\sigma}{\vec{\tau}}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}}$ is the electron spin operator. $\epsilon_{n}$ in the last line are the single particle eigenvalues for the electron system in a spin background $\{{\vec{m}}_{i}\}$, and $\mu=U/2$.

Calculating ${\partial F_{SF}}/{\partial\vec{m}_{i}}$, which is $\langle{\vec{s}}_{i}\rangle$, requires knowledge of the eigenvalues and eigenfunctions of the whole system. At $U/t_{hop}=3$ we can calculate $\langle{\vec{s}}_{i}\rangle$ accurately by constructing a cluster around the site ${\bf R}_{i}$ and diagonalising the cluster Hamiltonian instead of having to diagonalise the full $H$. We use a 13 site cluster centred on ${\bf R}_{i}$, including 4 nearest neighbour sites and 8 next nearest neighbour sites. We have benchmarked this against the full diagonalisation based calculation.

The electron temperature profile obtained from MFD can be approximated as: $T_{el}(t)=T_{el}^{i}e^{-t/\tau_{el}}+T_{el}^{f}(1-e^{-t/\tau_{el}})$. We find $\tau_{el}\sim 2$ ps and $T_{el}^{f}\propto T_{el}^{i}$ and a ratio $T_{el}^{i}/T_{el}^{f}=2.5$ allows us to model all $E_{0}$ using only one $\gamma$ value (see Supplement C). This leave $T_{el}^{f}$ as the only fluence dependent parameter in the Langevin equation. The fit $T_{el}^{f}(E_{0})$ is shown in the Supplement, where we also estimate $\gamma\approx 0.1$. We set bath temperature $T_{b}=40$K, the temperature for experiments were $\sim 80-100$K. We use the Euler-Maruyama algorithm to solve the LD equation with step size $\delta t\sim 0.01t_{hop}^{-1}$. We average the Langevin data over 5-10 runs when extracting timescales for the magnetic response.

Timescales: In Fig.1 we show the time dependence of the system averaged magnetic moment $m(t)=\frac{1}{N}\sum_{i}|\vec{m}_{i}|$ and the 2D structure factor $S_{\bf Q}(t)=\frac{1}{N^{2}}\sum_{ij}e^{i{\bf Q}.({\bf R}_{i}-{\bf R}_{j})}{\vec{m}}_{i}(t)\cdot{\vec{m}}_{j}(t)$ at ${\bf Q}=(\pi,\pi)$. We set the pre-pulse values of $S_{\bf Q}$ and $m$ to 1. As the pulse hits the system both $m(t)$ and $S_{\bf Q}(t)$ are suppressed on a timescale $\sim 0.3$ps. Post pulse, the mean moment rises quickly (panel (a)) while $S_{\bf Q}$ (panel (b)) has a strongly fluence dependent recovery time. The time axis in panel (b) is logarithmic, highlighting the wide range of recovery times. We fit the suppression-recovery dynamics in $m(t)$ and $S_{\bf Q}(t)$ to simple exponentials of the form: $I(t)=I(0)e^{-t/\tau_{d}}+I(\infty)(1-e^{-t/\tau_{r}})$ For both $m$ and $S_{\bf Q}$ we find $\tau_{d}\sim 0.2$ ps. Panel (c) shows the recovery time $\tau_{m}$ of $m(t)$, while panel (d) shows the recovery time $\tau_{ord}$ for $S_{\bf Q}$. $\tau_{m}$ ranges from $1-3$ps while $\tau_{ord}$ ranges from $1-120$ps. Size dependence of $\tau_{ord}$ is shown in Supplement D.

The associated long term amplitudes $m(\infty)$ and $S_{\bf Q}(\infty)$ are shown in panel (e). For an experimental system in a thermal environment, $m(t)$ and $S_{\bf Q}(t)$ should return to their pre-pump value after the post pump system attains equilibrium. Within our MFD this does not happen since a fraction of the pump created double occupancy persists at long time. Their deexcitation requires multimagnon emission processes which have a long timescale and are beyond MFD. Even in the experiments some indicators do not return to pre-pump values over experimental observation time.

Optical response: The most dramatic effect, and readily measurable consequence, of the excited electron population is in the optical conductivity $\sigma(\omega)$. At equilibrium the Mott insulator at low temperature should have no weight in the real part of $\sigma(\omega)$, for $\omega<\Delta$, the equilibrium gap. Since the electronic timescale is $\sim$ ten times shorter than magnetic timescale (as inferred from the magnetic bandwidth in Fig.3, later), we calculate $\sigma(\omega,t)$ using the instantaneous electronic eigenstates and eigenvalues at time $t$ (see Supplement E).

Fig.2 shows features of $\sigma(\omega,t)$ at $E_{0}^{2}=0.12$, the fluence at which $S_{\bf Q}$ first drops to zero, Fig.1.(b). In Fig.2(a) $\sigma_{R}(\omega)$ at $t=0$ shows the absence of any weight for $\omega<\Delta$. Between $t=0$ and $1$ps the weight in the interval $\omega=[0,\Delta]$ rises quickly, reaching a maximum around 1ps and declining thereafter. $\sigma(\omega)$ deviates from a Drude response due to the suppressed low energy density of states and the strong orientational disorder in the magnetic background.

In panel (b) we show the integrated weight over a small window $\omega=[0.05\Delta-0.1\Delta]$, as had been done in the pump probe experiment [30]. The two noteworthy features are (i) the quick decay to a long time value with a time constant $\tau_{opt}\sim 2$ps, (ii) a ‘long time’ value that is $\sim 50\%$ of the peak. Panel (c) shows the time dependence of the weight extracted from the experimental data in Fig.2(g) in [30].

It is not coincidental that the optical timescale $\tau_{opt}$ and the moment recovery time $\tau_{m}$ are comparable. Both arise from the deexcitation of electrons pushed to the upper band by the pump pulse. The deexcitation reduces double occupancy, trying to restore the moment magnitude, and the reducing number of electrons in the UHB, and ‘holes’ in the LHB, reduces the low energy optical weight. Broadly speaking, these process are related to the quick - local - charge relaxation in the system, operative on a few ps timescale. This contrasts with the long timescale for restoring global order. We next look at an indicator - already probed experimentally - of the momentum and ‘time resolved’ magnetic fluctuation spectrum.

Spin dynamics: Experiments have measured the 2D dynamical structure factor $S({\bf q},\omega)$ of the spins via time resolved resonant inelastic X-ray scattering (tr-RIXS) [29]. In contrast to the equilibrium case this collects ‘${\bf q}-\omega$’ data of magnetic fluctuations over a time interval $\pm\Delta t$ around a reference time $t$. The background state is time evolving so the resulting $S({\bf q},\omega)$ depends on the reference time $t$. We compute the corresponding object based on LD data: $S({\bf q},\omega,t)={1\over N^{2}}{\Big{|}}\sum_{i}\int_{t-\Delta t}^{t+\Delta t}dt^{\prime}e^{i{\bf q}.{\bf R}_{i}}e^{i\omega t^{\prime}}{\vec{m}}_{i}(t^{\prime}){\Big{|}}^{2}$ We set $\Delta t=10$ps and used $t=10,30,50,70$ps.

Fig.3(a)-(b) show $S({\bf q},\omega)$ as colour maps for two values of reference time, $10$ps and $70$ps. On the $x$ axis is ${\bf q}$, on the $y$ axis is $\omega$, and the intensity is coded in colour. The results are at $E_{0}^{2}=0.4$, a ‘strong fluence’ case where magnetic order recovers very slowly. Superposed on (a) and (b) is the equilibrium spin wave dispersion $\omega_{0}({\bf q})$ of the $(\pi,\pi)$ AF state. At both $t$ the spectrum far from the ordering wavevector $(\pi,\pi)$ looks similar to $\omega_{0}({\bf q})$ except for a correction due to moment size. However near $(\pi,\pi)$ the spectral weight distribution is very different from $\omega_{0}({\bf q})$: at $t=10$ps the difference is drastic while at $t=70$ps $S({\bf q},\omega)$ still has not recovered its character. Note that these times are much greater than $\tau_{m}$ or $\tau_{opt}$ that we have discussed.

To highlight the ${\bf q}$ dependence of the recovery process panels (c) and (d) show the lineshapes for two momenta: one situated far from $(\pi,\pi)$ at $(0,\pi)$, the other near $(\pi,\pi)$ at $(0.9\pi,0.9\pi)$, In each of these we have highlighted $\omega=\omega_{0}({\bf q})$ by an arrow. It is obvious that for $(0,\pi)$ the peak location in $S({\bf q},\omega)$ has stabilised within $10$ps. At $(0.9\pi,0.9\pi)$ however the peak location is far from stabilised even at $70$ps. This is not surprising given that the ordering timescale $\tau_{ord}$ at this fluence is $\sim 60-70$ps, Fig.1(d), and it would take $\sim 2-3$  $\tau_{ord}$ for the spectrum near ${\bf Q}$ to stabilise.

More generally, the results in Fig.3 suggest that the ‘recovery’ time for the magnetic fluctuation spectrum is strongly ${\bf q}$ dependent, and this timescale $\tau_{fluc}({\bf q})$ ranges from $\tau_{m}$ when ${\bf q}$ is far from ${\bf Q}$ to $\tau_{ord}$ when ${\bf q}\rightarrow{\bf Q}$. The RIXS spectrum probes moment size recovery and short range correlations when far from ${\bf Q}$, and progressively longer range correlation of the moments when ${\bf q}\rightarrow{\bf Q}$. In the next figure we wish to show the ‘domain growth’ physics that underlies this phenomena.

Spatial behaviour: Fig.4 shows the detailed spatial behaviour of the moment size $\vec{m}_{i}$ in the top row, differently oriented AF domains (coded by colour) in the middle row, and the instantaneous structure factor $S_{\bf q}$, in the bottom row, for time ranging from $-1$ps to $100$ps. This is a strong pulse situation, the same as studied for $S({\bf q},\omega)$.

As the pulse hits the system the mean magnetic moment reduces to 40$\%$ of its original value and then quickly recovers to a stable value by 10ps. There is no spatial structure to the fluctuations in $\vec{m}_{i}$. The middle row shows the pre-pulse perfect order at $t=-1$ps - a single AF domain. At $3$ps, where the moment value is still suppressed, there is a patchwork of small domains with linear dimension of a few lattice spacings (the map is $60\times 60$). By $10$ ps the moments have stabilised and there are only a few large competing domains. $50$ps and $100$ps show the increasing dominance of one (green-blue) domain. The lowest row shows the ${\bf q}$ dependence of $S_{\bf q}$ near the ordering wavevector. The perfect order in the pre-pulse state is destroyed at 3ps due to suppression of the moments. Then there is a slow growth of intensity near ${\bf q}={\bf Q}$ and a clear peak becomes visible only at $50$ps.

In the domain pictures we have gauged out the $(\pi,\pi)$ oscillation of Neel order and plotted equivalent ‘ferromagnetic’ domains with net moment pointing in different directions. If a domain has linear dimension $L_{d}$, probing it should lead to a fluctuation spectrum mimicking the bulk order as long as $q_{x},q_{y}\ll\pi/L_{d}$. The typical $L_{d}$ for us is time dependent, so the moderate size domains at $10$ps well capture the fluctuation at $(\pi,0)$, while we would need $t\gtrsim 100$ps to capture fluctuations near $(\pi,\pi)$.

We discuss two issues now that have bearing on the theory-experiment comparison and the reliability of the calculation itself. (i) 2D versus 3D: The experimental paper [29] suggested that the growth of 3D order was slow because the interplanar magnetic exchange $J_{\perp}$ was $\ll J_{||}$, the in plane value. It was implicit that recovery of 2D order was quick and the delayed 3D recovery was a $J_{\perp}/J_{||}$ effect. We have shown that 2D recovery itself is inevitably delayed. Does that mean an even more delayed 3D recovery? We did dynamics on a ‘layered’ Heisenberg model (Supplement F) with $J_{\perp}/J_{||}$ down to $10^{-3}$. Naively one may have expected $\tau_{3D}/\tau_{2D}\sim 10^{3}$. We find the ratio to be 2! This is because 2D is the lower critical dimension for $O(3)$ models and any small $J_{\perp}$ is a singular perturbation [43, 44]. (ii) Thermalisation: Within the MFD framework, from which we extract our $T_{el}(t)$, there are pump induced holon-doublon excitations which persist to long time, and their effective temperature differs from the apparent temperature sensed by the magnetic moments. For a system that equilibriates, these two temperatures should finally be the same. The electronic excitation scale $\Delta\gg J_{||}$ the magnetic excitation scale, so multimagnon emission processes are needed to deexcite electrons. The timescale arising from such processes has been estimated to be $\sim t_{hop}^{-1}e^{\alpha(\Delta/J_{||}log(U/t_{hop}))}$, [45] where $\alpha\sim O(1)$. Plugging in $\Delta\sim 500$ meV and $J_{||}\sim 25$ meV, we would get a decay time $10^{4}$ ps, which is $\gg$ than the experimental observation time or our run time.

Conclusion: Pump-probe experiments have made it possible to temporarily ‘metallise’ an antiferromagnetic Mott insulator and suppress it’s magnetic order by creating double occupancy and destroying the magnetic moment. The key question was how this strongly perturbed state evolves back towards the reference AFMI at long times. With the experimental data on Sr₂IrO₄ setting a reference we set up a hierarchical scheme that addresses this nonequilibrium correlated problem on large spatial scales in real time. We conclude that the ‘charge physics’, of optics etc, is dominantly local, quick, and mostly insensitive to fluence. The magnetic order recovery, however, is strongly non local, involves growth of domains, and brings in a fluence and system size dependent timescale. Momentum resolved magnetic excitations probe different spatial scales, and hence different recovery times. While our specific results are on the AFMI in Sr₂IrO₄, the mean field dynamics framework, and it’s reduced Langevin counterpart, can be readily adapted to address large spatial scale nonequilibrium phenomena in ordered systems like superconductors or charge density waves.

We acknowledge use of the HPC clusters at HRI. PM thanks Rajdeep Sensarma for a discussion.