Witnessing light-driven entanglement using time-resolved resonant inelastic X-ray scattering

Characterizing and controlling entanglement in quantum materials is crucial for the development of next-generation quantum technologies. However, defining a quantifiable figure of merit for entanglement in macroscopic solids is theoretically and experimentally challenging. At equilibrium the presence of entanglement can be diagnosed by extracting entanglement witnesses from spectroscopic observables and a nonequilibrium extension of this method could lead to the discovery of novel dynamical phenomena. Here, we propose a systematic approach to quantify the time-dependent quantum Fisher information and entanglement depth of transient states of quantum materials with time-resolved resonant inelastic x-ray scattering. Using a quarter-filled extended Hubbard model as an example, we benchmark the efficiency of this approach and predict a light-enhanced many-body entanglement due to the proximity to a phase boundary. Our work sets the stage for experimentally witnessing and controlling entanglement in light-driven quantum materials via ultrafast spectroscopic measurements.


SUPPLEMENTARY NOTE 1: DERIVATION OF THE SELF-CONSISTENT EQUATION
The nonequilibrium dynamic spin structure factor S(q, ω, t) defined in Eq. (2) of the main text reads as where we use the Wigner transformation of the time variablesτ = (t 1 +t 2 )/2 and τ = t 1 −t 2 .
In this paper, we employ a pure initial state |ψ(−∞) at zero temperature for all simulations, with the aim to reduce the computational complexity and focus on the nonequilibrium aspect of the problem. However, the expressions for the QFI and trRIXS spectra, as well as the relation between them [i.e. Eq. (6) of the main text] are not restricted to a pure state, but can be generalized to any mixed initial state with time-independent distribution weights.
TheÛ (t 1 , t 2 ) is the unitary time evolution operator, Z is the partition function of the equilibrium state (at t = −∞), and β is the inverse temperature. To simplify the derivation, we define the momentum-space spin excitation operator This equality states that the energy integral of the nonequilibrium dynamic spin structure factor at a given momentum q is a convolution of the two-time spin correlation functions with the time-dependent profile of a probe with finite pulse duration. Note that the QFI f Q (q, t) is defined as 4 ρ s −q (t)ρ s q (t) /N , when the SU(2) symmetry is preserved. (The disconnected part of Eq. (1) in the main text is non-zero in the presence of long-range magnetic order; however, it can be evaluated separately from the elastic scattering peak intensity and subtracted off from the correlation function). The expression for f Q (q, t) is an equal-time measurement for an instantaneous wavefunction at time t [1][2][3]. Therefore, without time-translational invariance in nonequilibrium systems, one cannot obtain the correct time evolution of the QFI by just integrating the dynamic spin structure factor.
However, the information about the time sequence of the QFI {f Q (q, t)| − ∞ < t < ∞} is encoded in the entire time evolution of S(q, ω, t), although without snapshot-to-snapshot correspondence. Therefore, one can apply self-consistent iterations to deconvolve f Q from the integral equation. By introducing a change of variable x = τ − t we can Taylor expand Noting that the integral on the right-hand side is nonvanishing only for even values of m, we can rewrite this equation in the form of Eq. (5) of the main text where Since simulating correlated electrons using the Krylov-subspace requires very small time steps, S(q, ω, t) is evaluated over a fine time grid and enables the reliable calculation of high-order derivatives.
Supplementary Eq. (6) can be numerically solved by a self-consistent iteration scheme where we truncate the derivative series to some finite m = M and start with the snapshot QFI as Then the k-step iteration depends on the lower orders as The convergence is reached by satisfying the criterion |f In practice, we employ the Fourier method to solve this iterative equation, which does not require an artificial termination.

HOLE LIFETIME
Due to the finite core-hole lifetime, RIXS does not reflect the exact S(q, ω). As shown in the Fig. 5 of the main text, trRIXS captures the oscillation of QFI for the wavefunction dynamics, but there is a finite offset between the values extracted from trRIXS and those evaluated exactly using the wavefunctions. It is important to note that this offset does not increase in time compared with that at equilibrium. This phenomenon indicates that finite lifetime effect, which causes the deviation between RIXS and S(q, ω), is insensitive to whether the initial state is equilibrium.
Therefore, we introduce a correction of the trRIXS intensity by a constant factor determined by comparing the equilibrium RIXS spectrum and the S(q, ω) (which can be measured with inelastic neutron scattering). As shown in Supplementary Fig. 1 for τ core = 1/2t h , the QFI extracted from equilibrium RIXS spectrum aligns with that from the S(q, ω) after scaling by an overall factor 1.163. After this rescaling, the QFI extracted from trRIXS is sufficiently accurate throughout the entire dynamics. The efficiency of this correction in turn demonstrates that trRIXS probes the spin excitations out of equilibrium with limited errors. This statement holds not only for the excitation energy, but also for the spectral intensity after accounting for this constant correction factor.

LUTION OF THE LOCAL MOMENT
To further test the self-consistent approach, we apply Eq.
As shown in Supplementary Fig. 3, the equilibrium local moment at quarter filling is n ∼  Supplementary Fig. 2, the compression of spin fluctuations at large momenta, inherent from the equilibrium Luttinger instability, is partially compensated by the enhancement at small ones (q = π/6 in this system) in the EHM. Such an enhancement is minor for the Hubbard model, as shown in the main text Fig. 6c. Therefore, the presence of the nonlocal interactions and the proximity to a phase boundary are crucial for the light-induced entanglement.

TROPY
In this section, we compare the nonequilibrium QFI dynamics with the entanglement entropy, which has been widely employed as an entanglement measure in quantum information. We consider the 1D extended-Hubbard model with the same parameters as the main text. The Rényi entropy for an instantaneous nonequilibrium state is defined as Here we consider the case α = 1. The label A refers to the label of the subsystem A, defined as half of the chain, and A c refers to the remaining subsystem.
In order to compare the time-dependent entanglement entropy with the QFI dynamics, we adopt the same pump condition as Fig. 4 of the main text, i.e. A 0 = 1 and Ω = 10t h . As shown in Supplementary Fig. 4, the nonequilibrium entanglement entropy increases at the same time of the light-driven QFI, further supporting the notion that the nonequilibrium state exhibits light-enhanced entanglement. It is, however, worth noting that, unlike the QFI, the entanglement entropy is not experimentally accessible in spectroscopic measurements of quantum materials.

SINGLE-PARTICLE FERMIONIC MODES
In this appendix, we follow Ref. 4 and derive the QFI bounds based on single-particle fermionic modes (SPFM). For the 1D extended Hubbard model of N sites (with periodic boundary conditions) considered in the main text, the operator ρ s q in Supplementary Eq. (3) obtained as a linear superposition of local operators S z i is non-hermitian at all momenta (ρ s † q = ρ s −q except q = 0, π). The corresponding Hermitian operator can be chosen as where s ↑ = +, s ↓ = −. The selection of local density operators in O q avoids ambiguities related to fermion anticommutation and sign. The fluctuation of this O q operator is related to the (instantaneous) QFI of the spin operator ρ s q through For the simulated translational-symmetric system in Sec. IV, the time-dependent wavefunction |ψ(t) has conserved momentum quantum number. This implies that (ρ s q ) 2 = 0 only when q = 0, π. At the same wavevectors q = 0, π, ρ s q is Hermitian, therefore Supplementary Eq. (12) leads to To estimate the upper bound on f Q (q, t), we consider the upper bound on O 2 q as derived in Ref. 4 for a fermionic many-body k-producible pure quantum state. For the creation (annihilation) operators c † nσ (c nσ ), all SPFMs are spanned by n = 1, 2, · · · N and spin σ =↑, ↓. If P m = [(n m 1 , σ m 1 ), (n m 2 , σ m 2 ), · · · , (n m Nm , σ m Nm )] denotes a unique partition that contains N m SPFMs, then P = P 1 ∪ P 2 ∪ · · · ∪ P Nparts (where N parts ≤ 2N is the number of partitions) forms the complete set of SPFMs. With this notation, Ref. 4 shows that if |Ψ(t) k−prod has a fixed number of electrons N e , then the lower bound on O 2 q = k−prod Ψ(t)|O 2 q |Ψ(t) k−prod can be determined by choosing a partitionP =P 1 ∪P 2 ∪ · · ·P Nparts which maximizes the right-hand side of the following inequality where A β (q) ∈ {σ cos(qn)|σ = ±1, n = 1, 2, · · · , N } such that A 1 (q) ≥ A 2 (q) ≥ · · · ≥ A 2N (q) andP m = P high m ∪P mid m ∪P low m represents a decomposition of the partitionP m . The right-hand side in Supplementary Eq. (14) can be numerically evaluated via the algorithm presented in Ref. 4, assuming conserved particle number N e . Thus, the upper bound of QFI based on