Nonadiabatic Exchange-Correlation Potential for Strongly Correlated Materials in the Weak and Strong Interaction Limits

: In this work, nonadiabatic exchange-correlation (XC) potentials for time-dependent density-functional theory (TDDFT) for strongly correlated materials are derived in the limits of strong and weak correlations. After summarizing some essentials of the available dynamical mean-ﬁeld theory (DMFT) XC potentials valid for these systems, we present details of the Sham–Schluter equation approach that we use to obtain, in principle, an exact XC potential from a many-body theory solution for the nonequilibrium electron self-energy. We derive the XC potentials for the one-band Hubbard model in the limits of weak and strong on-site Coulomb repulsion. To test the accuracy of the obtained potentials, we compare the TDDFT results obtained with these potentials with the corresponding nonequilibrium DMFT solution for the one-band Hubbard model and ﬁnd that the agreement between the solutions is rather good. We also discuss possible directions to obtain a universal XC potential that would be appropriate for the case of intermediate interaction strengths, i.e., a nonadiabatic potential that can be used to perform TDDFT analysis of nonequilibrium phenomena, such as transport and other ultrafast properties of materials with any strength of electron correlation at any value in the applied perturbing ﬁeld.


Introduction
Description of the physical properties of strongly correlated materials, i.e., systems with partially filled, localized d and f electron orbitals, is one of the greatest challenges of modern condensed matter physics and material science. A strong competition between the role of the usually itinerant s-and/or p-states, and the much more localized d-and/or f-states, in these systems produces many new effects that already have applications in technologies. The properties can be further tuned by changing pressure, temperature, the size and composition of the system and other parameters [1][2][3][4]. The nonequilibrium properties of these materials are of a special interest since they can be used in modern technologies that are pushing the limits of "small and fast". In particular, the ultrafast properties of strongly correlated materials can be used in smart window applications [5], biosensors [6] and switches [7]. It is thus imperative that reliable and accurate strategies be developed for both understanding and predicting the equilibrium and nonequilibrium properties of strongly correlated materials. This task is, however, challenging because of complex processes in these systems driven by both kinetic and potential energy (local electron repulsion) effects: technically, one cannot use a perturbation expansion in either of the corresponding (kinetic or potential energy) terms of the Hamiltonian, and rather, nontrivial, nonperturbative methods have to be applied.
There are also technical difficulties in applying the available nonperturbative manybody tools, such as equilibrium [8] and nonequilibrium [9,10] dynamical mean-field theory (DMFT), to real materials since they typically contain a good number of nonequivalent atoms and require the inclusion of several orbitals in the calculations and/or necessity GGA potentials, the reason could also be a large reduced gradient in the strongly correlated cases (of localized charges), which leads to noncorrect predictions [38]).
The situation is similar in higher dimensions. In particular, Karlsson et al. [47] obtained with DMFT the static dependence E XC (n) and the TDDFT adiabatic XC potentia V XC (n) = ∂E XC ∂n = V XC (n) → V XC [n](r, t) for the 3D one-band Hubbard model. Similar to the 1D case, it was found that the obtained XC potential does not work well when U is large and/or when the system is close to half-filling, once more suggesting importance of nonadiabatic effects (for the importance of these effects in the cases of Hubbard dimer and Hubbard chains, see, e.g., Refs. [56,57,59] and [60], respectively). Recently, the linearresponse TDDFT nonadiabatic XC potential (XC kernel) was calculated with DMFT for both spin-independent (infinite-dimensional Hubbard model) [61] and spin-dependent (noncollinear) (bulk Ni) [62] cases. In these works, the kernels were derived from the DMFT charge and spin susceptibilities, respectively (it must be mentioned that almost two decades ago, nonadiabatic XC kernels for Hubbard clusters were also obtained from charge susceptibility [54]).
To summarize, for "large" (bulk) systems, currently we have available nonlinear adiabatic (one-band Hubbard model) and linear nonadiabatic (many-band Hubbard model) DMFT XC potentials (more accurately, in the last case, tools to calculate them), but the problem of finding general TDDFT nonlinear nonadiabatic XC potentials for multi-band systems still needs to be solved (the situation is schematically shown in Figure 1, in which the connections between different approximations for the XC potentials are also shown). model, became inaccurate or even wrong as U increased or the doping approached halffilling (i.e., when the system approaches the regime of strong correlations) [33]. One of the reasons for this could be lack of nonadiabatic effects, i.e., missing memory in the XC potential (in the case of GGA potentials, the reason could also be a large reduced gradient in the strongly correlated cases (of localized charges) , which leads to noncorrect predictions [38]).
The situation is similar in higher dimensions. In particular, Karlsson et al. [47] obtained with DMFT the static dependence E (n) and the TDDFT adiabatic XC potential V (n) = = V (n) → V [n](r, t) for the 3D one-band Hubbard model. Similar to the 1D case, it was found that the obtained XC potential does not work well when U is large and/or when the system is close to half-filling, once more suggesting importance of nonadiabatic effects (for the importance of these effects in the cases of Hubbard dimer and Hubbard chains, see, e.g., Refs. [56,57,59] and [60], respectively). Recently, the linear-response TDDFT nonadiabatic XC potential (XC kernel) was calculated with DMFT for both spin-independent (infinite-dimensional Hubbard model) [61] and spin-dependent (noncollinear) (bulk Ni) [62] cases. In these works, the kernels were derived from the DMFT charge and spin susceptibilities, respectively (it must be mentioned that almost two decades ago, nonadiabatic XC kernels for Hubbard clusters were also obtained from charge susceptibility [54]). To summarize, for "large" (bulk) systems, currently we have available nonlinear adiabatic (one-band Hubbard model) and linear nonadiabatic (many-band Hubbard model) DMFT XC potentials (more accurately, in the last case, tools to calculate them), but the problem of finding general TDDFT nonlinear nonadiabatic XC potentials for multi-band systems still needs to be solved (the situation is schematically shown in Figure 1, in which the connections between different approximations for the XC potentials are also shown).   In this work, we derive nonadiabatic XC potentials by solving the Hubbard model exactly in the limits of weak and strong correlations when the system is perturbed by an external pulse. The derivation is based on an alternative tool-the Sham-Schluter equation approach [63,64], where V XC is obtained from the many-body theory results for the electron self-energy Σ ij t i , t j . We approximate the obtained self-energies by the local-in-space siteindependent functions, Σ ij t i , t j = δ ij Σ t i , t j , relevant to DMFT. In the many-body case, such an approach with nonequilibrium perturbative self-energy was successfully applied to study, e.g., the nonequilibrium current [65]. The obtained potential allows one to accurately take into account the correlation and memory effects in the cases of small and large Us. Additionally, the results can be regarded as a step towards the final goal of deriving the DMFT nonadiabatic XC potential for strongly correlated materials.
The rest of the text is organized as follows. In Section 2, we mention key facts about TDDFT and summarize the previous main results for the DMFT XC potential for the linear-response case for spin-nonpolarized and spin-polarized systems (Sections 2.1 and 2.2, correspondingly). In Section 3, we summarize the previous results for the nonlinear adiabatic XC potential obtained with DMFT. Details of the Sham-Schluter equation formalism are given in Section 4. In Sections 5 and 6, we present our results for the XC potentials obtained at small and large values of the local Coulomb repulsion U by solving the Sham-Schluter equation. We conclude, with a summary and a discussion of possible future directions in the field, in Section 7.

Linear Response: DMFT Nonadiabatic XC Kernels
To analyze the dynamics of a system with TDDFT, one needs to solve the timedependent Kohn-Sham (KS) equation [13,14]: where α are quantum numbers, except spin σ, and is the KS potential, a unique functional of density that consists of the Hartree V H [n](t, → r ), the XC V XC [n](t, → r ) and the external V(t, → r ) (ions, external fields, etc.) parts. Since one or more parts of the KS potential in Equation (1) depend on the charge density, this equation has to be solved self-consistently with the charge density equation: ( In the linear-response regime, the equations are simplified since the dependence of the XC potential on the change density n(t, → r ) ≡ n 0 ( → r ( + δn(t, → r ) (basically, the only unknown time-dependent quantity in the problem) is rather simple: where is the XC kernel. This quantity can be found from the TDDFT equation for the charge susceptibility:  (6) i.e., (in frequency representation), provided the susceptibility χ(t → r , t → r ) is known. The remaining undefined function in Equation (7) is KS susceptibility: where ε j and ϕ 0 j → r are the KS eigenenergies (band spectrum) and eigenfunctions, and In the case of solids, i.e., "infinite" periodic systems, it is more convenient to solve Equation (6) or (7) in the frequency-momentum representation. Namely, in the frequency representation one gets from Equation (6): Due to periodicity of the response function in real space, ( → R is the lattice vector), one can use the Fourier transforms for all functions of type where → G are reciprocal vectors (and similar for other functions). Then, Equations (6) and (7) in the frequency-momentum space can be written as: where v→ is the Fourier transform of the Coulomb potential, and are the matrix elements of the KS susceptibility. Provided the XC kernel is known, one can study the dynamics of the system by using Equations (1) and (3) and to analyze other physical properties of the system, such as the absorption spectrum, by using results for the TDDFT susceptibility in Equation (12) (see, e.g., [14]). Below, we discuss solutions of the opposite problem: finding the XC kernel from the results for susceptibility.

Nonadiabatic DMFT XC Kernel: Paramagnetic Case
In this subsection, we show how to derive the nonadiabatic XC kernel f XC (tr, t r ) from the Hubbard model susceptibility obtained in the DMFT approximation in the paramagnetic case [61]. For this, we begin by writing down the multi-orbital Hubbard Hamiltonian in standard notation: where t lm ij are the inter-site, inter-orbital hopping parameters, µ is the chemical potential, J is the exchange energy parameter and n l iσ are the site-and orbital-dependent spin density. For the system described by this Hamiltonian, one can calculate the charge susceptibility for Equation (6) by using the following equation: where n(r, t) is the total charge density operator andT is the time-ordering operator. The functions in the last sum are two-particle spin-orbital susceptibilities: The noninteracting ("KS") susceptibility χ 0 for Equation (6) (and Equation (13)) can also be obtained from Equation (17), but at U = J = 0.
In DMFT, it is assumed that both one-electron GFs and susceptibilities are local in space. We do not present here the details on how to solve the DMFT problem but just mention that DMFT is based on the approximation of momentumindependent electron self-energy, i.e., local-in-space self-energy that corresponds to an effective one-site problem and is exact in the limit of infinite dimensions, or infinite atomic coordination number. Instead, we give details on the main steps of the calculation of susceptibility and, hence, of the XC kernel. These details are instructive since they also form the basis of the many-body calculation of the XC kernel by using approaches other than DMFT.
To calculate the one-frequency susceptibility, one can use the Bethe-Salpeter equation approach, beginning with the calculation of the generalized four-time local susceptibility (in the imaginary time representation, path integral form; for details, see Ref. [61]): where G is the dynamical mean field, a function that describes the effects of the rest of the system on the electrons on the given site, in the effective one-site problem (the last time argument in the correlators can be set to zero due to the time-translation invariance).
After the most complicated step-calculation of the susceptibility (20)-is finished, one can obtain its Fourier transform: and substitute the result into the equation for the local vertex function Γ(ω m , ω l , ω n ): (we skip indices for simplicity), where are two-frequency susceptibilities. Then, from Equations (22) and (23), one can calculate the vertex functions Γ(ω m , ω l , ω n ) and use the result to calculate the needed DMFT susceptibility: χ(ω n ) is used to calculate the DMFT XC kernel from Equation (13), or, more explicitly: (we use the real frequency representation in the last equation). One can also use a simplified, semi-analytical, approach to calculate the susceptibility with the ladder approximation for the vertex function [66][67][68]: where U ch is an effective (charge-channel) interaction parameter. The resulting formula for the susceptibility is: where is familiar one-loop susceptibility. To find the unknown parameter U ch , one can solve the equation that connects the susceptibility with the average charge density and the double occupancy: Computation 2022, 10, 77 is the one-site electron double occupancy (for the multi-orbital case, see Ref. [68]). The ladder formalism was applied to the one-band half-filled Hubbard model with [61]. It was found that the kernel demonstrates oscillations with frequencies on the order of typical excitation energies (between Hubbard sub-bands and states involving zero-energy quasi-particle peak). The kernel also demonstrates the memory time of the order of the inverse hopping parameter. An example of the time dependence of the DMFT XC kernel is shown in the bottom-right of Figure 1 (for details, see Ref. [61]). Importantly, a similar structure of the local kernel was found in the cases of 2D 3 × 3 Hubbard cluster and Hubbard dimer in Ref. [54] (in this work, the results are in frequency representation).

Nonadiabatic DMFT XC Kernel: The Spin-Polarized Case
The approach presented in the last subsection was also generalized on the spindependent case [62], including the case of noncollinear systems with nonconserved total spin (e.g., systems with helical spin-wave ground states and systems with varying surface magnetization). In the noncollinear spin (spin-flip) TDDFT approach, the KS equation has the following form: where the wave functions depend on spin and the potentials are spin matrices. In particular, in this equation V XCσσ [n](r, t). is a functional of the spin-density matrix In the linear-response case, where is the XC kernel matrix. Similar to the paramagnetic case, in the DMFT approximation the non-Hartree part of where f DMFT XCσσ σσ (t − t ) is the DMFT local XC kernel matrix that satisfies the equation In the last equation, we introduced the local susceptibility matrix (in time representation) Computation 2022, 10, 77 9 of 23 (and the corresponding noninteracting susceptibility matrix χ (0) σσ σσ (t)) and expressed all matrices in a 4 × 4 matrix representation with independent indices for spin-index pairs 1 =↑↑ , 2 =↑↓ , 3 =↓↑ and 4 =↓↓ . In the one-loop approximation for the susceptibility (37) it was found that the XC kernel matrix has the following structure: where the nondiagonal matrix elements are responsible for the spin-flipping processes [62].
The frequency-and time-dependencies of the matrix in the Equation (38) elements are semi-quantitatively rather similar to the dependencies of the XC kernel in the paramagnetic case, including~1 fs memory time (bottom right Figure 1). The XC kernel matrix (38) was used to analyze the ultrafast spin dynamics and the ultrafast demagnetization in Ni after perturbation by a laser pulse [62]. It was shown, in particular, that spin-flip processes generated by the correlation effects taken into account in the DMFT XC kernel can lead to an ultrafast demagnetization in the system, in agreement with experimental data (see Ref. [62]).

Nonlinear Response: Adiabatic XC Potential
Karlsson et al. [47] used DMFT with the Lanczos exact diagonalization algorithm for the impurity problem to calculate the XC energy and the XC potential for the 3D oneband Hubbard-model Hamiltonian (16), with an additional term, ∑ i,σ ε i n iσ , to include the possible effects of nonhomogeneity (with the energy ε i and also the Coulomb repulsion U i that depend on the site). From the DMFT solution for the one-particle GF, the authors constructed the energy functional [52]: In Equation (39) and for the potential energies were used to calculate the DMFT ground state energy of the system Due to some noise in the solution for the XC energy at low doping and at close to half-filling, a smoothed polynomial fitting at 0.2 < n < 1 was used to obtain the analytical result [40] at low doping: (λ is a fitting parameter).
DMFT XC potential V XC [n] was obtained by differentiating E XC [n] with respect to n. Due to particle-hole symmetry, the potential is antisymmetric with respect to point n = 1, i.e., V XC [2 − n] = −V XC [n]. It was found that when U reaches the value of the MIT, U = 14t, the XC energy curve E XC [n] has a cusp at half-filling that gives a discontinuity in the XC potential (in the 1D case, the discontinuity in the XC potential (obtained with the 1D case Bethe ansatz LDA potential [20,21]) exists at any U, i.e., there is no MIT). Typical doping dependencies of the obtained DMFT XC potential in the metallic and insulating phases are shown in the cartoon top-left of Figure 1.
As it was mentioned in Section 1, the derived XC potential was used to study response of a system (a 125-atom cubic system with a non-zero Coulomb interaction at the central atom, a system that can also be solved exactly) to a pulse perturbation, and it was found that the agreement between the exact results and the DMFT-TDDFT results differ in the cases of strong correlations at large Us and/or close to half-filling. Since the potential is adiabatic, it was suggested that memory effects might be important in these regimes. Another reason for the difference in the results might be not being included in the XC potential nonlocal effects. Thus, while the local DMFT approximation is rather good in the static case, one cannot exclude that nonlocal effects might be important when the system is out of equilibrium. Indeed, recently it was demonstrated [25] that a nonlocal XC self-energy correction to the DMFT XC potential (obtained with a time-perturbative many-body approach) significantly improves the DMFT XC (adiabatic LDA) potential results.
Similarly, the Hubbard-Holstein model for the electron-phonon system was solved with DMFT, and the DMFT XC energy and the XC potential were calculated [48]. It was found that in the infinite-dimensional case at large Us, as in the Hubbard model, the XC potential has a cusp at n = 1. It was also found that the interaction between electrons and phonons delays the inset of the discontinuity as the interaction U increases and that this interaction also reduces the height of the cusp.

Sham-Schluter Equation Formalism
The Sham-Schluter equation approach [63,64] is a method that allows one to obtain, in principle, exact nonadiabatic and nonlocal XC potentials from many-body theory. Before presenting our results for the XC potentials obtained with this approach, we present the main details of this formalism, first in the one-band case and then in the multi-band case.
In the one band case, one can write down the Dyson equation for the electron GF G and the self-energy Σ: where G 0 is the free-electron GF that satisfies the equation (v(t 1 → r 1 ) is the external (ion, perturbation, . . . ) potential).
In the differential form, Equation (45) reads (in this equation, the time arguments are defined on the complex Kadanoff-Baym time contour [10], and the time-dependencies of different quantities are obtained by using the time arguments on the first "physical time" branch, also used in TDDFT).
An equation similar to Equation (47) can be written for the KS GF: . (48) (with an equation similar to Equation (46) for the free-electron KS GF). Combining Equations (47) and (48), one can obtain an equation that connects the many-body self-energy and V XC : where is the XC self-energy. Since the many-body and the KS GF give the same charge density, one can obtain from Equation (49): This is the Sham-Schluter equation that connects the (TD)DFT XC potential V XC and the many-body self-energy Σ. It can be transformed to a more convenient form: where (it is important that Equation (53) is not the solution for V XC (t 1 → r 1 ), since χ and G KS in this equation are functionals of V XC ).
The system of Equations (47), (48) and (52) that one needs to solve to find the XC potential is very complicated. Thus, some approximations have to be made. One of them could be G(t 1 → r 1 , t 2 → r 2 ) ≈ G KS (t 1 → r 1 , t 2 → r 2 ), where the time-ordered KS GF is: In Equation (55), φ l (t 1 → r 1 ) (l defines the state number) are the wave functions obtained from the solution of the KS equation. Then, the problem reduces to solving the KS equation for φ l (t 1 → r 1 ) with the XC potential defined by Equations (52)- (54). The solution gives both φ l (t 1 → r 1 ) and the XC potential (an alternative approach, such as by Krieger-Li-Iafrate (KLI) [69], can be also used to solve this TDDFT version of the optimized effective potential problem [70]).
The methodology for the multi-orbital is rather similar, with the many-body orbital GF G lm (t 1 → r 1 , t 2 → r 2 ) equation used instead of Equation (47).
Here again, one repeats the steps for the one-orbital case to obtain: Using again n(t 1 , , one can obtain the Sham-Schluter equation generalized for the multi-orbital case [61]: where Thus, to find V XC from Equation (59), one needs two functions: the KS GF and the many-body multi-orbital self-energy (the remaining functions G lm can be found from the self-energy by using Equation (56)). Yet, in this case, the problem is even more complicated than in the one-band problem and approximations have to be made. Again, one of them is: where the KS GF G KS (t 1 → r 1 , t 2 → r 2 ) is defined in Equation (55) (and can be obtained by solving the KS equation). The self-energy that defines the XC potential can be found from the solution of a many-body problem. In particular, in the case of DMFT, this function has the following structure: Naturally, it is difficult to solve the nonequilibrium many-body problem for Σ lm (t 1 → r 1 , t 2 → r 2 ). Therefore, it is desirable to have a simple analytical result for this function. Below, we derive analytical formulas for the nonequilibrium self-energy for the one-band Hubbard model in the cases of small (Section 5) and large (Section 6) values of U.

Nonadiabatic XC Potential: Small-U Solution
To obtain the solution for the second order in U self-energy for the Hubbard model, we calculate the contour time-ordered GF by using the perturbation-theory approach: whereĤ pert is the perturbation part of the HamiltonianĤ pert . In the case of small Us, H pert (t) = U ∑ in ↑i (t)n ↓i (t). To obtain connection between the self-energy and GF in different orders of approximation, let us write the Dyson equation: and expand the GF and the self-energy in this equation in different orders of the perturbation theory: (66) (in this Section, the expansion is actually performed in powers of dimensionless parameter U/t, and one needs to restore the hopping parameter where it is appropriate). From the last equation one can obtain the equations that connect the first-and the second-order self-energy and the corresponding GFs: Thus, obtained from Equation (64) G (1) and G (2) allows us to obtain Σ (1) and Σ (2) from Equations (67) and (68). After rather simple calculations for the nonequilibrium self-energy in the second order in U, one can obtain: The corresponding formula in the momentum domain is (we write down the spin-up expression, but the result for the spin-down case can be easily obtained from this equation) ↓k +q (t 2 , t 1 )G (0) ↓k+q (t 1 , t 2 ). (70) The last result gives the following formula for the equilibrium self-energy in the frequency representation (after the Fourier transformation with respect to t 1 − t 2 ): The self-energy in Equation (71) corresponds to the vertex-corrected Σ(ν) that gives momentum-dependent correction to this function in DMFT [71]: where, in our case, the vertex is F ↑↓kk q (ν, ν , ω) = −U. The momentum-dependent correction to the DMFT self-energy is an important hot topic in modern research. In this work, we are focusing on the local GFs and self-energies, in the spirit of DMFT. In this case, The next step is to substitute the corresponding matrix GFs and self-energies, instead of the contour time-ordered functions in Equation (73). In Equation (74), are the retarded, advantage and Keldysh GFs, respectively (these functions are needed to calculate the corresponding time-ordered functions on the first branch of the contour, see below). From the matrix form of Equation (73), one obtains: To be explicit, we consider the case of the hypercubic lattice relevant to DMFT with the spectrum ε(k) = −2t ∑ i cos(k i ) and the corresponding Gaussian free-electron DOS t * 2 , where t * is the rescaled hopping (we also consider the paramagnetic case). In the case of an external perturbation, the hopping matrix elements are modified by the Peierls substitution: We assume that the field is homogeneous in space and the vector potential points in the direction of the diagonal of the unit cell. In this case, the free-electron spectrum becomes: (we drop the factor 1 c in front of the vector potential) with the two-energy DOS ρ 2 (ε, ε) = 1 πt * 2 a d e with and Σ > XC (t 2 , t 3 ) and Σ < XC (t 2 , t 3 ) are defined in Equations (98) and (99), respectively. Thus, we have derived the expression for the orbital XC potential in the case of small Us. The following further simplifications can be made: (i) Space-local approximation (in spirit of DMFT), where all space vectors under the integral are the same, and (ii) to dispose of the pulse dependence of the self-energy (to make the potential "perturbation-invariant"), one can neglect the correlation effects during the pulse and put in Equations (98) and (99) S(t 1 , t 2 ) = 0, C(t 1 , t 2 ) = t 1 − t 2 (these results follow from Equations (88) and (89) at A(t) = 0). This gives semi-analytical results: The exponential decay with time of the XC self-energies in the above lends itself to extraction of the memory time (e.g., for the Gaussian DOS it can be found performing the energy integration in the last two equations and taking into account the factors e − 1 To give a feeling of the structure of the XC self-energy and hence of the XC potential, in Figure 2a we show results for Σ > (t 1 , t 2 = 0) for several pulse durations. We find that, for short pulses, the structure of Σ > (t 1 , t 2 = 0) is quite similar after the pulse time (Figure 2a), suggesting that the approximation A(t) = 0 for the XC potential might be meaningful.

Nonadiabatic XC Potential: Large-U Solution
In the case of large U, perturbation theory, Equations (64)-(68), can be also applied, but with H (t) = − ∑ t c (t)c (t) , , (the expansion is performed in powers of dimensionless parameter t/U). Here, the main difference from the small U case is that the perturbation part of the Hamiltonian includes time explicitly, as modified by the Peierls substitution (82) for the hopping matrix elements.
In the first-order approximation, one then obtains (we consider the case with spindiagonal GFs and self-energies):

Nonadiabatic XC Potential: Large-U Solution
In the case of large U, perturbation theory, Equations (64)-(68), can be also applied, but withĤ pert (t) = − ∑ l,m,σ t lm c + lσ (t)c mσ (t) (the expansion is performed in powers of dimensionless parameter t/U). Here, the main difference from the small U case is that the perturbation part of the Hamiltonian includes time explicitly, as modified by the Peierls substitution (82) for the hopping matrix elements.
In the first-order approximation, one then obtains (we consider the case with spindiagonal GFs and self-energies): and Σ (1) Since we assume that the atomic GFs G ilσ (t, t 1 ) are local in space, one obtains: and Σ (1) i.e., the linear in-hopping part of the self-energy is zero.
One thus obtains: (we restored the dependence on U). The atomic GF does not depend on hopping, thus there is no vector potential in this function and, hence, there is no time-dependence. Therefore, G Next, in the matrix representation, Equation (114) gives: where G R(0) ii,σ (t 1 , Computation 2022, 10, 77 18 of 23 and the other two GFs are: From Equations (115)-(120), one can obtain: and Together with Equations (98) and (99) for the small U case, Equations (124) and (125) for the case of strong interaction are the main results of this article. These equations define the XC potential (100) in the case of large Us. As in the case of weak interactions, here one can also calculate a spatial-local and a pulse-independent approximation, which gives cos(eA(t 1 ) − eA(t 2 )) = 1. However, contrary to the case of small U, the XC potential (124), (125) demonstrates oscillating behavior, with no signs of a finite memory time. The oscillations can be eliminated by using the adiabatic approximation, Σ(t 1 , t 2 ) ∼ δ(t 1 − t 2 ). However, this approximation gives a rather trivial XC potential (100) with time-independent time-ordered self-energy.
In Figure 2b, we show the time dependence of the greater self-energy (125) in the cases of different pulse durations. Again, as one can see from this figure, the details of the pulse become nonimportant at times longer than pulse, suggesting that the approximation A(t) = 0 is a meaningful one in this case as well.
To test the kernel, we applied it to the half-filled one-band 3D Hubbard model at U = 4t * ≡ 4 √ 6t (in the more challenging case of a large U) in presence of an ultrafast laser pulse and compared the result for the time-dependent excited charge density with the corresponding results obtained with nonequilibrium DMFT (nonequilibrium iterativeperturbation theory solver). As it follows from the comparison displayed in Figure 3, the agreement between the results of nonequilibrium DMFT and DMFT-TDDFT is rather good. This is not surprising, since in both theories the time-dependence of the interacting dynamics is defined by basically the same function-nonequilibrium electron self-energy.
corresponding results obtained with nonequilibrium DMFT (nonequilibrium it perturbation theory solver). As it follows from the comparison displayed in Figu agreement between the results of nonequilibrium DMFT and DMFT-TDDFT i good. This is not surprising, since in both theories the time-dependence of the int dynamics is defined by basically the same function-nonequilibrium electron self

Summary and Future Directions
In this work, we have derived an analytical form of the time-ordered self-en the Hubbard model for both weak (Equations (95), (98) and (99)) and strong (Eq (95), (124) and (125)) local Coulomb repulsion. The obtained results were used to c a nonadiabatic XC potential for strongly correlated systems (with an approximat given in Equation (100)) by using the Sham-Schluter equation approach (below, notations V and V for the weak and strong interaction potentials, respective only assumption used in the derivation of the self-energies was the locality of the the spirit of DMFT. Now that we have obtained potentials that correspond to the cases of w strong Coulomb interactions, the logical question is: how do we deal with the int case of intermediate U? The first attempt could be to use a hybrid potential V = (1 − x)V , where 0 ≤ x ≤ 1. The value of the parameter x would depend on whe system is closer to the weak (x ≪ 1) or strong (x~1) correlation regimes or is b them (x~0.5) (a more rigorous approach to define a weak-strong regime can be lated in order to derive the value of the mixing parameter). A more rigorous app to solve the nonequilibrium DMFT problem and to find the self-energy in Equat and to use it to construct the XC potential in Equation (100). Moreover, such a n librium numerical solution of the DMFT problem to construct the XC potential b √ 6t when the system is perturbed by pulse (magenta curve). In the DMFT-TDDFT case, the atomic hydrogen s-orbital wave functions were used for the DFT input.

Summary and Future Directions
In this work, we have derived an analytical form of the time-ordered self-energy for the Hubbard model for both weak (Equations (95), (98) and (99)) and strong (Equations (95), (124) and (125)) local Coulomb repulsion. The obtained results were used to construct a nonadiabatic XC potential for strongly correlated systems (with an approximated form given in Equation (100)) by using the Sham-Schluter equation approach (below, we use notations V WI XC and V SI XC for the weak and strong interaction potentials, respectively). The only assumption used in the derivation of the self-energies was the locality of the GFs, in the spirit of DMFT. Now that we have obtained potentials that correspond to the cases of weak and strong Coulomb interactions, the logical question is: how do we deal with the interesting case of intermediate U? The first attempt could be to use a hybrid potential V XC = xV WI XC + (1 − x)V SI XC , where 0 ≤ x ≤ 1. The value of the parameter x would depend on whether the system is closer to the weak (x 1) or strong ( x ∼ 1) correlation regimes or is between them ( x ∼ 0.5) (a more rigorous approach to define a weak-strong regime can be formulated in order to derive the value of the mixing parameter). A more rigorous approach is to solve the nonequilibrium DMFT problem and to find the self-energy in Equation (95) and to use it to construct the XC potential in Equation (100). Moreover, such a nonequilibrium numerical solution of the DMFT problem to construct the XC potential by using Equation (100), or more generally Equation (59), is necessary in the case of lattices different from the hypercubic one and/or in the multiband case, where it is not possible to obtain a simple analytical expression for the self-energy. The exact sum rules for the nonequilibrium self-energy [72,73] can be used to test the accuracy of the obtained potential. Since the multi-orbital nonequilibrium DMFT problem is computationally intensive, an effective one-orbital model (with the same values of U, average hopping and charge density per orbital) might be used to obtain Σ XC (t 1 , t 2 ). We tested the potential for the one-band Hubbard model and the obtained results that are encouraging.
It is important to mention that the potentials obtained are orbital potentials. There is thus the need to use a large enough number of orbitals in each particular case, depending on the width of the frequencies of the applied pulse. A convenient way to solve the corresponding TDDFT problem is the density matrix formalism (see, e.g., Ref. [74]). On the other hand, due to an analytical form of the time-dependent functions in the integral in the XC potential (Equation (100)), the computational cost of the approach is relatively low, comparable to other TDDFT methods (in particular, the calculation of the result shown in Figure 3 took approximately 2 h).
Another important point is the sum rules. Though we use a quasi-local local-inspace approximation (more accurately, our potential is orbital, not LDA), it is known that even local (TDLDA) XC potentials satisfy a very important exact sum rule-the harmonic potential theorem [75]. On the other hand, it is also known that in the case of exactly local nonadiabatic potentials, TDDFT might give contradictory results [76]. Since the potential in Equation (100) is not local in space, this issue should not arise. Nevertheless, many more tests of this functional are needed to obtain a better idea about its accuracy.
We expect the potentials can be used to analyze different properties of strongly correlated systems, including for novel materials such as nanoparticles and twisted bilayer graphene, from optical response to superconductivity (for the current status in the field, see Ref. [77]; though, in the superconducting case in addition to the derived "normal diagonal" self-energies one needs to add "nondiagonal" (Gorkov) self-energy parts, and the XC potential should also have more components (as discussed, e.g., by Gross and collaborators in Ref. [78] in the static case).
Several important tasks still need to be accomplished: (1) To understand the relations of the obtained potentials with the adiabatic Karlsson-Privitera-Verdozzi DMFT XC potential [47] and with the nonadiabatic XC kernel [61] (and its extension to the spin-polarized case [62]) (see Figure 1). (2) To perform further analysis in order to understand the dependence of the XC potential on the external perturbation field and to make the potential "universal", i.e., fieldindependent (putting A = 0 in the expressions for the self-energy can be regarded as the first attempt acceptable in the case of short pulses). (3) To integrate the role of phonons in the results for the spin-polarized case (see, e.g., Ref. [48]) to be able to apply the method to longer times. (4) To generalize the approach to the nonlocal case, i.e., to take into account the momentumdependence of the XC self-energy [71].

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
All data that support the findings of this study are included within the article.