{'fraction_non_alphanumeric': 0.04445108770638826, 'fraction_numerical': 0.014061996964540312, 'mean_word_length': 4.576264666281977, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-time time-dependent density functional theory using higher order finite element methods

We present a computationally efficient approach to solve the time-dependent Kohn-Sham equations in real-time using higher-order finite-element spatial discretization, applicable to both pseudopotential and all-electron calculations. To this end, we develop an a priori mesh adaption technique, based on the semi-discrete (discrete in space but continuous in time) error estimate on the time-dependent Kohn-Sham orbitals, to construct a close to optimal finite-element discretization. Subsequently, we obtain the full-discrete error estimate to guide our choice of the time-step. We employ spectral finite-elements along with Gauss-Legendre-Lobatto quadrature to render the overlap matrix diagonal, simplifying the inversion of the overlap matrix that features in the evaluation of the discrete time-evolution operator. We use the second-order Magnus operator as the time-evolution operator in all our calculations. Furthermore, the action of the discrete Magnus operator, expressed as exponential of a matrix, on the Kohn-Sham orbitals is obtained efficiently through an adaptive Lanczos iteration. We observe close to optimal rates of convergence of the dipole moment with respect to spatial and temporal discretization, for both pseudopotential and all-electron calculations. We demonstrate a staggering 100-fold reduction in the computational time afforded by higher-order finite-elements over linear finite-elements. We present comparative studies, in terms of accuracy and efficiency, of our approach against finite-difference based discretization for pseudopotential calculations, and demonstrate significant computational savings when compared to the finite-difference method. We also demonstrate the competence of higher-order finite-elements for all-electron benchmark systems. Lastly, we observe good parallel scalability of the proposed method on many hundreds of processors.


Introduction
Time-dependent density functional theory (TDDFT) extends the keys ideas of ground-state density functional theory (DFT) to electronic excitations and time-dependent processes. It relies on the Runge-Gross theorem [1] to establish, for a given initial state, a one-to-one correspondence between the time-dependent external potential and the time-dependent electronic density, thereby making the electronic density the fundamental variable to define other physical quantities. Subsequently, one invokes the Kohn-Sham ansatz [2] to reduce the many-electron time-dependent Schrödinger equation to a set of effective single electron equations, called the time-dependent Kohn-Sham (TDKS) equations. Although this entails, analogous to the ground-state case, the use of approximate exchange-correlation functionals, TDDFT offers a great balance of accuracy and computational efficiency to enable the study of a wide-range of time-dependent phenomenaoptical [3] and higher-order responses [4,5], electron transport [6,7], charge-transfer excitations [8,9], dynamics of chemical bonds [10], multi-photon ionization [11,12,13], to name a few.
Over the past two decades there has been a growing interest in developing faster and more accurate numerical methods for solving the TDKS equations. Broadly, these numerical methods can be classified into two categories, characterized by the strength of the light-matter interaction, namely, linear-response time-dependent density functional theory (LR-TDDFT) [14,15] and real-time time-dependent density functional theory (RT-TDDFT) [16,17,18]. The LR-TDDFT pertains to the case of weak interaction between the external field and the system, wherein the field induces a small perturbation from the ground-state. In such perturbative regime, one can compute the linear density response from the ground-state itself, which in turn can be used for the calculation of first-order response functions such as the absorption spectra. The RT-TDDFT, on the other hand, is a more generic framework which captures the electronic dynamics in real-time, thereby, allowing to handle both perturbative and non-perturbative regimes (e.g., harmonic generation, electron transport) in a unified manner. This involves propagating the TDKS equations in real-time without any restriction to the external field in terms of its frequency, shape or intensity.
Despite its generality in dealing with various time-dependent processes, there are two major challenges associated with RT-TDDFT. The first stems from the quality of the time-dependent exchange-correlation approximation used in the TDKS equations. The exact exchange-correlation functional is, in general, nonlocal in both space and time [19,20,21] and has an initial-state dependence [22]. However, the lack of insight into its time nonlocality and initial-state dependence has necessitated the use of the adiabatic approximation, wherein the exchange-correlation functional is defined in terms of the instantaneous electronic density. Although the applicability of the adiabatic approximation to various systems and materials properties are yet to be understood, they have shown remarkable agreement in estimating the transition frequencies [3], and, in most cases, is the underlying approximation in existing RT-TDDFT packages. As with most of the numerical implementations in RT-TDDFT, this work is restricted to the adiabatic approximation. The second challenge stems from the huge computational cost associated with the non-linear TDKS equations. Numerical simulations for large length-and time-scales are still computationally challenging, and warrant systematically improvable, accurate, efficient and scalable spatio-temporal discretizations.
Significant efforts have been made towards efficient RT-TDDFT numerical schemes as extensions to popular ground-state DFT packages, borrowing from their respective spatial discretization. These include planewave basis in QBox [23,24]; linear combination of atomic orbitals (LCAO) in Siesta [25,26] and GPAW [27]; Gaussian basis in NWChem [28,29]; and finite-difference based approaches in Octopus [30] and GPAW [31,32]. The planewave basis, owing to its completeness, provides systematic convergence, and affords an efficient treatment of the electrostatic interactions through fast Fourier transforms. However, they remain restricted to only periodic geometries and boundary conditions, thereby ill-equipped to describe systems with defects, and non-periodic systems like isolated molecules and nano-clusters. Additionally, the nonlocality of the basis greatly hinders its parallel scalability. Atomic-type orbitals, such as LCAO and Gaussian basis, owing to their atom-specific basis, are well-suited to describe molecules and nano-clusters for both pseudopotential as well as all-electron calculations. However, owing to the incompleteness of such basis, systematic convergence for all materials systems remains a concern. Furthermore, the nonlocality of the basis substantially limits parallel scalability of their numerical implementations. The finite difference discretization (FD) provides systematic convergence, can handle a broad range of boundary conditions, and exhibits improved parallel scalability in comparison to planewave and atomic-type orbital basis. However, incorporating adaptive spatial resolution in FD through a non-uniform grid remains non-trivial. Additionally, the lack of basis in FD makes an accurate treatment of singular potentials difficult, thereby restricting its applicability to only pseudopotential calculations. On the other hand, the finite-element basis [33,34], being a local-piecewise polynomial basis, offers several key advantages-it provides systematic convergence; is amenable to adaptive spatial resolution, and thereby suitable for both pseudopotential and all-electron calculations; exhibits excellent parallel scalability owing to the locality of the basis; and admits arbitrary geometries and boundary conditions. The efficacy of the finite-element basis in terms of its accuracy, efficiency, scalability and relative performance with other competing methods (e.g., planewaves, Gaussian basis, FD), have been thoroughly studied in the context of ground-state DFT, for both pseudopotential [35,36,37,38,39,40,41,42,43,44,45,46] and all-electron calculations [35,45,46,47,48,49,50,51,52,53]. A similarly comprehensive study on the efficacy of the finite-element basis for RT-TDDFT is, however, lacking. While two recent studies [54,55] demonstrate the accuracy of finite-elements for RT-TDDFT, they remain restricted to only linear and quadratic finite-elements. As has been known from prior studies in ground-state DFT [56,48,45], linear and quadratic finite-elements require a large number of basis functions (50, 000 − 500, 000 per atom for pseudopotential calculations) to achieve chemical accuracy, and hence, perform poorly in comparison to planewaves and other real-space based methods. However, this shortcoming of linear and quadratic finite-elements for ground-state DFT calculations has been shown to be alleviated by the use of higher-order finite-elements [56,57]. In particular, the works of Motamarri et. al. [58,45] have demonstrated a 100-to 1000-fold improved efficiency of the higher-order finite-elements for ground-state DFT calculations over linear finite-elements. In this work, we extend the use of higher-order finite-elements to RT-TDDFT calculations and demonstrate the resulting advantages.
The keys ideas in this work can be summarized as: (i) developing an a priori mesh-adaption, through a semi-discrete (discrete in space, continuous in time) error analysis of the time-dependent Kohn-Sham equations, and subsequently, obtaining a close to optimal finite-element discretization for the problem; (ii) use of spectral finite-elements in conjunction with Gauss-Legendre-Lobatto quadrature to render the overlap matrix diagonal, thereby simplifying the evaluation of the inverse of the overlap matrix that features in the discrete time-evolution operator; (iii) obtaining a suitable temporal discretization using a full-discrete error analysis, in the context of second-order Magnus time-evolution operator; and (iv) using an adaptive Lanczos iteration to efficiently compute the action of the Magnus propagator on the Kohn-Sham orbitals. The a priori mesh-adaption in this work draws from the ideas described in [59,60,58,45], wherein the optimal spatial discretization is obtained by minimizing the discretization error in the observable of importance, subject to fixed number of elements in the finite-element mesh. In particular, we proceed by obtaining a semi-discrete error estimate in the dipole moment of the system as a function of the characteristic mesh-size distribution, h(r), and the exact time-dependent electronic fields. Next, we minimize this estimate with respect to h(r), subject to fixed number of elements in the mesh, to obtain the optimal spatial discretization. The resulting expression for the optimal h(r) is defined in terms of the degree of the interpolating polynomials and higher-order derivatives of the exact time-dependent electronic fields. We remark that although the exact electronic fields are a priori unknown, the use of asymptotic behavior of the atomic wavefunctions [61] provides a practically useful strategy to coarse-grain the mesh away from the nuclei, thereby efficiently resolving the interstitial and vacuum regions in our numerical studies. We note that the above semi-discrete analysis is also useful in decoupling the spatial discretization from the choice of time-integration scheme. Having obtained the optimal spatial discretization, the suitable temporal discretization is obtained through a full-discrete error analysis, in the context of second-order Magnus time-evolution operator. This work, to the best of our knowledge, is the first work that guides the spatio-temporal discretization for the RT-TDDFT problem using error estimates. In contrast to the aforementioned a priori scheme, an a posteriori scheme offers an alternative approach by defining the optimal spatio-temporal discretization in terms of the discrete solutions, instead of the exact solutions. However, such an approach involves error analysis based on simultaneous spatial and temporal discretization, as opposed to the semi-discrete approach described above, and thus, has to be treated differently for different time integrators. While schemes for a posteriori error control, in the context of nonlinear Schrödinger equation, have been developed [62,63,64], they all employ Galerkin schemes for time-integration, which violate the unitary requirement of the time-evolution operator. An a posteriori spatio-temporal adaption scheme for RT-TDDFT based on unitary time-integrators (e.g., Crank-Nicolson, Magnus) is, to the best of our knowledge, lacking, and presents a future extension of this work.
We study the key numerical aspects of the proposed higher-order finite-element discretization for benchmark systems involving both nonlocal pseudopotential and all-electron calculations. To begin with, we study the numerical rates of convergence of higher-order finite-element discretization of the RT-TDDFT problem. We remark that although optimal rates of convergence have been demonstrated in the context of nonlinear Schrödinger equation [65,66,67,68,69], a similar study for RT-TDDFT is missing. In this work, we use two benchmark systems: (i) a pseudopotential calculation on methane molecule; and (ii) an all-electron calculation on lithium hydride molecule, to demonstrate the rates of convergence for linear, quadratic and fourth-order finite-elements. We employ the previously mentioned a priori mesh-adaption scheme to generate the meshes for all the benchmark systems. We observe numerical rates of convergence in the dipole moment close to the optimal rate of O(h p ), where p denotes the order of the finite-element. Next, we as-sess the computational advantage afforded by higher-order finite-elements over linear finite-element, using the same benchmark systems. We observe an extraordinary 100-fold speedup in terms of the total computational time for the fourth-order finite-element over linear finite-element, for calculations in the regime of chemical accuracy. We also compare the relative performance of the finite-element discretization against finite-difference method for pseudopotential calculations, for other material systems. For the pseudopotential case, we use aluminum clusters (Al 2 and Al 13 ), and Buckministerfullerene (C60) molecule as the benchmark systems. The finite-difference based calculations are done using the Octopus package [30]. Based on the benchmark system, the finite-element discretization shows a 3 to 60-fold savings in computational time as compared to the finite-difference approach, for pseudopotential calculations. We also demonstrate the efficacy of finite-elements for systems subjected to strong perturbation by studying higher harmonic generation in Mg 2 . Additionally, we demonstrate the competence of finite-elements for all-electron calculations on two benchmark systems-methane and benzene molecule. Lastly, we study the strong scaling of our implementation and observe good parallel scalability with ∼ 75% efficiency at 768 processors for a benchmark system of a Buckministerfullerene molecule containing 3.5 million degrees of freedom.
The rest of the paper is organized as follows. In Section 2, we briefly discuss the time-dependent Kohn-Sham equations and the form of the exact time-evolution operator. In Section 3, we derive the semi-and full-discrete error estimates, and use them for the spatio-temporal discretization. In Section 4, we describe the various numerical implementation aspects pertaining to spectral finite-elements and the discrete second-order Magnus operator. Section 5 details the convergence, efficiency and parallel scalability of the higher-order finite-elements along with its relative performance against the finite-difference. Finally, we summarize our findings and outline the future scope in Section 6.

Time-dependent Kohn-Sham Equations
TDDFT relies on the Runge-Gross theorem [1] and the Kohn-Sham ansatz [2] to reduce the manyelectron time-dependent Schrödinger equation to a set of effective single electron equations, called the timedependent Kohn-Sham (TDKS) equations. These equations prescribe the evolution of an auxiliary system of non-interacting electrons that yield the same time-dependent electronic charge density, ρ(r, t), as that of the interacting system. The TDKS equations, in atomic units, is given as where H KS [ρ](r, t; R), V KS [ρ](r, t; R) and ψ α (r, t) represent the time-dependent Kohn-Sham Hamiltonian, potential and orbitals, respectively, with the index α spanning over all the N e electrons in the system. R = {R 1 , R 2 , ..., R Na } denotes the collective representation for the positions of the N a atoms in the system. The electron density, ρ(r, t), is given in terms of the Kohn-Sham orbitals as In the present work, we restrict ourselves to only non-periodic (clusters and molecules) as well as spinunpolarized systems. However, we note that all the ideas discussed subsequently can be generalized to periodic (solids) as well as spin-polarized systems.
The time-dependent Kohn-Sham potential, V KS [ρ](r, t; R) in Eq. 1, is given by where V ext (r, t; R) denotes the external potential, V H [ρ](r, t) denotes the Hartree potential, and V XC [ρ](r, t) represents the exchange-correlation potential. The exchange-correlation potential, V XC [ρ](r, t), in general, is nonlocal in both space and time [19,20,21], and has a dependence on the initial many-electron wavefunction [22]. However, in absence of the knowledge of its true form, most of the existing approximations use locality in time (adiabatic exchange-correlation) and non-dependence on the initial many-electron wavefunction. This allows for direct use of the existing exchange-correlation approximations used in ground-state DFT. In the present work, we use the adiabatic local-density approximation (ALDA) [70], which is local in both space and time. Specifically, we use the Ceperley-Alder form [71].
In Eq. 3, the Hartree potential is given by The external potential comprises of the nuclear potential V N (r; R) and the external field V f ield (r, t). The nuclear potential is given by where Z I and R I represent the atomic charge and position of the I th nucleus. For a typical pseudopotential calculation, V psp N comprises of a local part, V loc psp , and a nonlocal part, V nl psp . For the nonlocal part, the action on a function φ(r), written in the Kleinman-Bylander form [72], is given by where l and m denote the angular and magnetic quantum number, respectively. Y lm (r, R I ) denotes the spherical harmonics at r with the origin at R I , δV I l is the specified l angular component short-ranged potential for the atom at R I , and L I is the maximum angular quantum number specified for the atom at R I . The external field, V f ield (r, t), is typically provided as a monochromatic laser pulse of the form where E 0 represents the electric field and g(t) represents the time-dependent envelope.
We note that both the electrostatic potentials-Hartree and nuclear (all-electron)-are extended in real space. However, using the fact that the 1 |r| kernel in these extended interactions is the Green's function of the Laplace operator, one can recast their evaluation as the solutions to the following Poisson equations: In the above equation, b(r; R) = − Na I=1 Z I δ(r; R I ), where δ(r; R I ) is a bounded regularization of the Diracdelta distribution with compact support in a small ball around R I and satisfies δ(r; R I ) dr = 1. We refer to previous works on finite-elements based ground-state DFT calculations [37,39,42,45,58] for a detailed treatment of the local reformulation of the electrostatic potentials into Poisson equations.
Formally, the solution to Eq. 1 can be written as where U (T, t 0 ) represents the time-evolution operator (propagator) and T denotes the time-ordering operator. Although the above equation provides a formal way to directly evaluate the orbitals at any time, t, resolving the implicit time-dependence of the Kohn-Sham Hamiltonian through the density is too difficult. However, one can exploit the following composition property of the propagator, to accurately resolve the implicit time-dependence in H KS [ρ](r, t). To elaborate, the above property allows us to rewrite the propagator U (T, t 0 ) as where t N = T and t i+1 − t i = ∆t i , with ∆t i denoting the variable time step. Consequently, one can divide the evaluation of the orbitals at T into N short-time propagation, given by In addition to resolving the implicit time-dependence in H KS [ρ](r, t), the short time propagation provides the numerical advantage of containing the norm of the exponent in Eq. 9. To elaborate, any efficient numerical scheme to compute the action of the propagator on a wavefunction involves either a power series expansion or a subspace projection of the propagator, wherein the number of terms in the power series or the dimension of the subspace required for a given accuracy are dependent on norm of the exponent. Moreover, there is a physical upper bound imposed on the time step based on the maximum frequency, ω max , that one wants to resolve in their calculations, i.e., ∆t max = 1 ωmax . Typically, ω max is determined by the eigenspectrum of the ground-state Hamiltonian or by the frequency of the applied field, V f ield . We note that, in practice, one uses a time step ∆t ∆t max owing to the need of containing time-discretization errors that arise in approximating the continuous propagator, T exp −i t+∆t t H KS [ρ](r, τ )dτ . We discuss these approximations and their associated time-discretization errors in greater detail in Sec. 3.3.

Error estimates and a priori mesh adaptation
In this section we provide an a priori mesh adaption based on minimizing the discretization error in the Kohn-Sham orbitals. To this end, we decompose the discretization error in the Kohn-Sham orbitals into two parts, one arising due to spatial discretization and the other due to temporal discretization. To elaborate, if ψ h α (r, t) and ψ h,n α (r, t) represent the semi-discrete (discrete in space but continuous in time) and full-discrete solution to ψ α (r, t), respectively, then one can use the triangle inequality to bound the discretization error in ψ α as follows where ||.|| denotes an appropriately chosen norm.

Notations, assumptions and preliminaries
Before deriving formal bounds on the two right-hand terms in the above equation, we introduce some notations, key assumptions and preliminaries. For a bounded closed domain Ω and bounded time interval [0, T ], we denote Ω T = Ω × [0, T ]. For any two complex-valued functions f (r, t), g(r, t) : Ω T → C, the inner product (f, g)(t) = Ω f (r, t)ḡ(r, t) dr, whereḡ(r, t) denotes the complex conjugate of g(r, t). Correspondingly, the norm ||f || L 2 (Ω) (t) = (f, f )(t). Thus, we extend the definition of the standard L 2 (Ω) and H 1 (Ω) spaces to define where Df denotes the spatial partial derivatives of f , and ∂Ω denotes the boundary of Ω. Additionally, we define two more spaces relevant to the Poisson problem (Eq. 8), For conciseness of notation, in all our subsequent discussion, we drop the argument t from the inner product as well as all the L p and H 1 norms. Thus, any occurrence of (. , .), ||.|| L p (Ω) , and ||.|| H 1 (Ω) are to treated as time-dependent, unless otherwise specified.
We now list certain weak assumptions that we invoke throughout our error-estimates.

A1
The time-dependent Kohn-Sham orbitals and their spatial derivatives are bounded and have a compact support on Ω, which, in turn, is a large but a bounded subset of R 3 . To elaborate, A2 The nuclear potential (in the all-electron case) due to the regularized charge distribution b(r; R) (defined in Eq. 8b) is bounded, i.e., V ae N ∈ L ∞ (R 3 ). A3 The local part of the pseudopotential is bounded, i.e., V loc psp ∈ L ∞ (R 3 ).
To begin with, we note, through straightforward use of Cauchy-Schwarz and Sobolev inequalities, that Furthermore, for the convolution integral of ρ and 1 |r| , denoted by |r| −1 * ρ = Ω ρ(x) 1 |r−x| dx, we have where we have used the Young's inequality along with the fact that |r| −1 ∈ L 2 (Ω).
We now attempt to bound , by decomposing V KS into its Hartree (V H ), nuclear (V N ), exchange-correlation (V XC ) and field (V f ield ) components, and bounding each of the components. For the Hartree potential (see Eq. 4), we have, for ∀v ∈ H 1 0 (Ω T ), In order to bound the two terms on the right of the above equation, we make the following observations. Firstly, by definition, Second, using A1, we note that ψ 2,α ∈ L ∞ (Ω T ). Further, as a consequence of Eq. 17a, V H [ρ Ψ1 ] ∈ L ∞ (Ω T ). Thus, using these arguments, we have, ∀v ∈ H 1 0 (Ω T ), For the first inequality in the above equation, we have used Eq. 17b in conjunction with Hölder inequality.
Next, for the exchange-correlation potential, we begin with the observation that for a given ρ X defined by the orbitals {χ 1 , χ 2 , . . . , χ Ne }, we have, Further, invoking the mean-value theorem, we observe for some χ α = λ α ψ 1,α +(1−λ α )ψ 2,α , with λ α ∈ [0, 1]. Thus, using the above relation, we have, ∀v ∈ H 1 0 (Ω T ), where we have used the Cauchy-Schwarz inequality in the second line, the Hölder inequality in the third line, and have invoked the boundedness assumption on V XC and V XC (assumption A5) in the last inequality.
For the nuclear potential in the all-electron case, we have, ∀v ∈ H 1 0 (Ω T ), where we have used Cauchy-Schwarz and Hölder inequalities in the first and second line, respectively. In the last line, we have invoked the boundedness assumption on V ae N (assumption A2).
For the pseudopotential case, we have, ∀v ∈ H 1 0 (Ω T ), where we have, once again, made use of the Cauchy-Schwarz and Hölder inequality along with the boundedness assumptions on V loc psp and δV I l (assumptions A3 and A4, respectively).
Lastly, invoking the boundedness assumption on V f ield (assumption A6), we get, ∀v ∈ H 1 0 (Ω T ), We now define the weak solution of the TDKS equation (Eq. 1) as follows: given an initial state ψ α (r, 0) ∈ Similarly, the weak solutions to the Poisson problems defined in Eq. 8 are defined to be

Semi-discrete error estimate
We denote X h,p ∈ H 1 (Ω) to be the finite-dimensional space of dimension n h , spanned by finite-element basis functions of order p. Further, we denote X h,p 0 = X h,p ∩ H 1 0 (Ω). We now define the semi-discrete solution, ψ h α (r, t), to Eq. 26 as follows: given an initial state where We now elaborate on the different terms appearing in the expression for represents the discrete counterpart of the weak solution described in Eq. 27a, and is obtained as the solution to Similarly, V h N (r; R), for the all-electron case, represents the discrete analog of the weak solution defined in Eq. 27b, and is given by For the pseudopotential case, V h N (r; R) is same as the continuous function V psp N (r; R).
We now introduce the concept of Ritz projection, P h , which will be used in subsequent error estimates. The Ritz projection P h : T ] is defined through the following Galerkin orthogonality condition, This allows us to use some standard finite-element error estimates (see [73]) in order to bound ψ − P h ψ L 2 (Ω) . In the following error-estimates we take a triangulation T h,p of p th order finite-elements covering the domain Ω [73].
Proposition 3.1. Assuming uniqueness and existence of the solution to Eqs. 26 and 28, we obtain the following bounds on the finite-element semi-discrete approximation error to the Kohn-Sham orbitals: where e denotes a finite-element of mesh size h e and cover Ω e in the triangulation T h,p , h min represents the smallest element in the triangulation T h,p , and |.| p,Ωe is the semi-norm in H p (Ω e ). The arguments s 1,α , s 2,α , and s 3 are defined as T ] in Eq. 26 and subtracting it from Eq. 28, we get and derive bounds on each of the terms. For simpler notation, we use u α = ψ α − P h ψ α and w α = (P h ψ α − ψ h α ). Thus, using ψ α − ψ h α = u α + w α , we rewrite Eq. 34 as Taking Noting that 1 2 and comparing the imaginary parts of Eq. 36, we have In the above equation, we note that (∇u α , ∇w α ) = 0, as a consequence of Eq. 31. Furthermore, (∇w α , ∇w α ) is real. Thus, Eq. 38 simplifies to, We now decompose V KS into its components to rewrite the second term on the right of the above equation as We note that the term V N ψ h α − V h N ψ h α , w α , on the right side of the above equation, is relevant only in the all-electron case (i.e., zero for the pseudopotential case as V N = V h N ). Combining the results from Eqs. 19,22,23,24, and 25, with v = w α , and using the fact that ψ h ∈ L ∞ (Ω), it is straightforward to show that Using the above result in Eq. 39, we obtain where we have used Eq. 16a in the second line to simplify the term involving ρ − ρ h L 1 (Ω) . Thus, summing the above equation over all index α, we have where in the second line we have split ψ α − ψ h α into u α and w α . Now, integrating the above equation, gives Noting that u α = ψ α − P h ψ α , ∂uα ∂t = ∂ψα ∂t − P h ∂ψα ∂t , and using the definitions of s 1,α , s 2,α , and s 3 (see Eq. 33), we can simplify the above equation as Invoking the Grönwall's inequality on the above equation yields Expressing Plugging the above relation into Eq. 46, we obtain We now introduce formal bounds on the terms involving , by using standard finite-element error estimates. To begin with, we use the definition of the Ritz projection along with Ceá's lemma [73] to note that For the ∂uα ∂t L 2 (Ω) term, we obtain, by taking the time-derivative of Eq. 31, We emphasize that, in arriving at the above equation, we choose v h that is time-independent. The above equation allows us to invoke the Ceá's lemma (see [73]) on ∂uα ∂t , to obtain Furthermore, using Eq. 1, we have, where we have used the definition of the |.| p+3 semi-norm in the third line and the boundedness assumptions (Ω) terms, we make use of standard finite-element error estimates for the Poisson problem to get Thus, using Eqs. 49, 52, and 53 in Eq. 48, we have Finally, expressing ψ α − ψ h α = w α + u α and combining the result of the above equation with Eq. 49, we obtain This concludes the proof of Eq. 32a.
In order to derive estimates for Squaring the above the inequality and summing over all the elements e, it is straightforward to obtain Thus, taking the square-root on the above inequality, we have Additionally, we once again invoke the Ceá's lemma [33] on u α = (ψ α − P h ψ α ) to note that Thus, using Eqs. 58 and 59, we get Further, using the inequality obtained in Eq. 54 in the above equation, yields Using the fact that h e /h min ≤ C for all the elements in T h,p , we can simplify the above equation as This concludes the proof for Eq. 32b.

Full discrete error estimate
We now derive the full-discrete error estimate in the context of the second-order Magnus time-evolution operator. To begin with, we represent the finite-element approximation to the Kohn-Sham orbitals and the Hartree potential as where N j : 1 ≤ j ≤ n h denote a basis for X h,p . Subsequently, the weak formulation in Eq. 28 can be written in the form of the following matrix vector equation where H jk and M jk denote the discrete Hamiltonian and overlap matrix, respectively, and ψ α (t) denotes the vector containing the coefficients ψ j α (t). The Hamiltonian matrix is, in turn, split into two parts, U jk and V jk , given by We note that the matrix M, owing to its positive definiteness, guarantees the existence of a unique positive definite symmetric square root, denoted by M 1/2 . This allows us to rewrite Eq. 67 as whereψ α (t) = M 1/2 ψ α (t) and {H,Ū,V} = M −1/2 {H, U, V}M −1/2 . In the following discussion, we drop the superscript h in ρ, for the sake of brevity.
We now invoke the Magnus ansatz to write the solution of Eq. 70 as The operator exp (A(t)) is termed as the Magnus propagator. Differentiating the above equation and com-paring with Eq. 70 gives [74]ψ with ad 1 X (Y) = XY − YX, and ad 0 X (Y) = Y. Now, inverting the operator on the right side of the above equation provides the following differential equation for A(t) [74]: where β k is given by the series e z −1 Finally, invoking the Picard iteration on the above equation yields the following explicit form of the Magnus expansion [75], From the above expression, we note that A is skew-Hermitian, owing to the Hermiticity ofH. As mentioned in Section 2, we resolve the implicit dependence ofH(t) on ρ(r, t) by using the composition property of a propagator (see Eq. 11). Thus, the Magnus propagator can be rewritten as where A n is given by Eq. 75 with the limits of integration modified to [t n−1 , t n ]. In practice, one replaces the exact A n with an approximate operatorÃ n , which involves, first, a truncation of the Magnus expansion, and second, an approximation for the time integrals in the truncated Magnus expansion. Truncating the Magnus expansion (defined in Eq. 75) after the first p terms results in a time-integration scheme of order 2p. In this work, for simplicity, we restrict ourselves to the second-order Magnus propagator, i.e., obtained by truncating the Magnus expansion after the first term.
We now express the time-discrete approximation toψ α (t n+1 ), denoted byψ n+1 α , as In particular, the action of the second-order Magnus propagator, with a mid-point integration rule, on a set of Kohn-Sham orbitals Φ = {φ 1 , φ 2 , . . . , φ Ne } which defines the density ρ Φ , is given by where ∆t = t n − t n−1 andH ρ t n−1 + ∆t 2 is the time-continuous Kohn-Sham Hamiltonian described by the Kohn-Sham orbitals Φ at the future time instance t n−1 + ∆t/2. We remark that the operator e −iH[ρ(tn−1+ ∆t 2 )]∆t results from two successive approximations: a) truncation of the Magnus expansion after the first term; and b) from the use of the mid-point integration to evaluate the time integral. Furthermore, H ρ t n−1 + ∆t 2 being dependent on a future instance of the density, is evaluated either by an extrapolation ofH using m(> 2) previous steps or by a second (or higher) order predictor-corrector scheme.
In the following error-estimates, we assume each time interval [t n−1 , t n ] to be of length ∆t. Moreover, for simple terminology, we term eÃ n as the second-order Magnus propagator without explicitly spelling out the mid-point integration rule invoked in it. Proposition 3.2. For a second-order Magnus propagator with a mid-point integration rule, we obtain the following bound for the time-discretization error to the vectorsψ α Proof. To begin with, we introduce the following operators To elaborate, S k 0 denotes the exact Magnus propagator from t 0 to t k , and R n k denotes the second-order Magnus propagator from t k to t n . We now rewrite the time-discretization error inψ α at t n using the following telescopic series, Noting that S k 0 = e A k S k−1 0 and R n k−1 = R n k eÃ k , we rewrite the above equation as Since R n k is a unitary operator, bounding ψ α (t n ) −ψ n α reduces to finding the bound on e A k − eÃ k ψ α (t k−1 ). To this end, we extend the proof presented in [75] to the non-linear case of the TDKS equations. We let a α k = e A k − eÃ k ψ α (t k−1 ) and split it as The two terms on the right hand side of the above equation denote the error due to truncation of the Magnus expansion and the time integral approximation, respectively.
In order to bound the error in the first term on the right side of Eq. 83, we introduce the following auxiliary function where B k (t) = t t k−1 −iH[ρ(τ )]dτ . We remark that ξ k α (t) encapsulates the effect of truncation of the Magnus expansion in the time interval [t k−1 , t k ]. Differentiating the above equation giveṡ whereḠ k (t) = idexp B k (t) (Ḃ k (t)). We observe thatḠ k is Hermitian. This can be proven as follows. First, note that for two Hermitian (or skew-Hermitian) matrices X, Y, the operator ad X (Y) is skew-Hermitian. Second, owing to the Hermiticity ofH, both ) are skew-Hermitian. Thus, by expanding dexp B k (t) (Ḃ k (t)) and using the above two arguments, it can be shown thatḠ k is Hermitian. Letting γ k α (t) =ψ α (t) − ξ k α (t), we use the Eqs. 72 and 85 to write the time-derivative of γ k α aṡ Now, taking the dot product with γ k α (t) † on both sides yields We note that 2Re{γ k where . represents the Euclidean norm of a vector. Further, is real, owing to the Hermiticity ofḠ. Thus, comparing the real parts of the above equation results in 1 2 Consequently, Time integrating the above equation yields where we have used the fact that γ k ) and the definition of the operator dexp X (Y), we get In order to bound [H(τ ),H(σ)], we begin by rewritingH in terms ofŪ andV, i.e., its kinetic and Kohn-Sham potential components. Noting thatŪ is time-independent, we Taylor expandH(σ) about τ to get whereV (τ ) = d dt (V(σ))| σ=τ . Thus using the above relation in Eq. 91 we get We now invoke the boundedness assumption onV (assumption A8), and the norm equivalence ofŪ andH (assumption A7), to obtain We now turn to the definition of ξ k α (defined in Eq. 84) to note that where the second equality follows from the fact that B k (t k ) = . We use this relation in Eq. 90 to get Plugging in the result of Eq. 94 into the above equation provides the following bound We note that the above equation provides a bound for the first term (truncation error) on the right side of Eq. 83. In order to bound the second term on the right side of Eq. 83, i.e., the error due to mid-point quadrature rule, we begin with the following identity Furthermore, we note that for a function f (x) if F 1/2 denotes the midpoint approximation to Using this result along with the unitarity of the operators e (1−x)Ã k and e xĀ k , we obtain , and invoking the boundedness assumption on d 2 dt 2V (assumption A8), we get Finally, using the results of Eqs. 97 and 100 in Eq. 82 along with the unitarity of the operators R n k , yields This concludes the proof of Eq. 79.

Optimal spatio-temporal discretization
We follow on the lines of [60,58,45] to obtain the optimal spatial discretization by minimizing the semidiscrete error in the dipole moment at a given time, subject to a fixed number of finite-elements. The dipole moment is a vector-valued function defined at time t by µ(t) = Ω x ρ(t) dx. We remark that the choice of dipole moment as an observable for this exercise is solely a matter of convenience, and any observable which can be inexpensively evaluated in terms of the density or the Kohn-Sham orbitals can be used instead.
Representing the x-component of the continuous and the semi-discrete dipole moment at time t as µ x (t) and µ h x (t), respectively, we use the Cauchy-Schwarz inequality to obtain where we have used the fact that x ∈ L 2 (Ω), for a finite domain Ω . Plugging Eq. 16b into the above equation, we have Now using Eq. 64b in the above equation, we have |ψ α | p+1,Ωe (s 1,α ) + |ψ α | p+1,Ωe (s 2,α ) + |ψ α | p+3,Ωe (s 2,α ) where s 1,α , s 2,α , and s 3 are defined in the Proposition 3.1. Using the definition of semi-norms, we rewrite the above equation as where N elem denotes the number of elements in the finite-element mesh (triangulation), and D p+1 and D p+3 denote the (p + 1)th and (p + 3)th spatial derivative of a function, respectively. Taking an element size distribution function, h(r), the above inequality can be further rewritten as, Obtaining the optimal element size distribution, for a fixed N elem , now reduces to the following minimization problem, Solving the Euler-Lagrange equation for the above objective function with respect to h(r) yields the following optimal element size distribution, where the constant E is evaluated from the constraint on the number of elements. We remark that although the above discretization approach requires a priori knowledge of the unknown ψ α (s 1,α ), ψ α (s 2,α ), and V H [ρ h ](s 3 ), we use the knowledge of their asymptotics from ground-state atomic calculations to coarsegrain the finite-element mesh for all calculations. Although, this does not inform us about the optimal discretization required near the nuclei, this affords an efficient strategy to accurately handle the regions away from the nuclei, wherein the time-dependent Kohn-Sham orbitals, typically, have similar decay properties.
Having determined the required spatial discretization, a suitable temporal discretization for the given finite-element mesh can be estimated by using the time-discrete error bound for the dipole moment. To elaborate, if µ h n,x denotes the x-component of the full-discrete dipole moment then using the result in Eq. 101 it is straightforward to show, As is evident from the above relation, our choice of ∆t is intrinsically tied to the spectral radius of the Kohn-Sham Hamiltonian, through ψ α Ū . Furthermore, we remark that although the presence of t n in the above inequality indicates increasing time-discretization error with time, it does not pose a limitation in containing the errors due to the fact that t n ≤ T , where, typically, T lies between 10 − 30 fs. Now, in order to evaluate a suitable ∆t that can contain the above error bound to a fixed tolerance, we need to estimate the values of ψ α (t) Ū and the value of the constant, C, featuring in it. The value of ψ α (t) Ū can be reliably approximated from its ground-state value, i.e., ψ α (0) Ū . The characteristic value of the constant is determined from atomic RT-TDDFT calculations at different ∆t. To elaborate, the constant can be evaluated from a linear fit to the log-log plot of the error |µ h x (t n ) − µ h n,x | with respect to ∆t. For a multi-atom system, we use the least ∆t obtained for each of its constituent atomic species.

Numerical Implementation
We now discuss some of the key numerical aspects involved in our implementation of the finite-element discretization of the TDKS equations.

Higher-order spectral finite-elements
Finite-elements, with their varied choices of forms and orders [76,34], have been widely used in several engineering applications. While the use of linear finite-elements remains popular in engineering applications that warrant moderate levels of accuracy, it remains computationally inefficient to attain chemical accuracy in electronic structure calculations. To highlight, the use of linear finite-elements have been shown to require large number of basis functions per atom (∼ 100, 000) to achieve chemical accuracy in ground-state DFT calculations [56,48]. However, this shortcoming has been demonstrably mitigated by the use of higher-order finite-elements [58,45]. In this work, we explore the possibility of similar gains from using higher-order finite-elements for RT-TDDFT calculations. Unlike conventional finite-elements, we employ spectral finiteelements for the spatial discretization of the TDKS equations. To elaborate, while the conventional finiteelements are constructed from a tensor product of Lagrange polynomials interpolated through equidistant nodal points in an element, spectral finite elements employ a distribution of nodes obtained from the roots of the derivative of Legendre polynomials or the Chebyshev polynomials [77]. In our work, we use the roots of the derivative of Legendre polynomials along with boundary nodes, so as to maintain C 0 continuity. The resulting distribution is known as the Gauss-Legendre-Lobatto node distribution. The spectral finiteelements are known to afford better conditioning with increasing polynomial degree [77], and have been used to gain significant computational efficiency in ground-state DFT calculations [58]. Moreover, a major advantage of the spectral finite-elements is realized when used in conjunction with the Gauss-Legendre-Lobatto (GLL) quadrature rule for evaluating the integrals involved in the overlap matrix M, wherein the quadrature points are coincident with the nodal points in the spectral finite-elements. Such a combination renders M diagonal, thereby greatly simplifying the evaluation of M −1/2 that features in the discrete TDKS equations (see Eq. 70). We note that while an n point rule in the conventional Gauss quadrature rule integrates polynomials exactly up to degree 2n − 1, an n point GLL quadrature rule integrates polynomials exactly only up to degree 2n − 3. Therefore, we employ the GLL quadrature rule only in the construction of M, while the more accurate Gauss quadrature rule is used for all other integrals featuring in the discrete TDKS equations. We refer to Motamarri et. al. [45] for a discussion on the accuracy and sufficiency of GLL quadrature in the evaluation of overlap matrix M. For the sake of brevity, we use the term finite-elements instead of spectral finite-elements in all subsequent discussions.

Approximating the second-order Magnus operator
The form of the Magnus operator, as shown in Eq. 77, calls for efficient means of evaluating exp(Ã n )ψ n α . Direct means of evaluating the matrix exp(Ã n ) remain computationally prohibitive beyond small sizes of A n . Alternatively, one can attempt to evaluate the action of exp(Ã n ) on ψ n α in an iterative fashion. Several such schemes are in use in RT-TDDFT calculations, namely, Taylor series expansion, Chebyshev polynomial expansion, split-operator techniques, and Krylov subspace projection method. We refer to Castro et. al. [78] and references therewithin for a detailed discussion on each of these schemes. In this work, we adopt the Krylov subspace projection method for its superior efficiency and robustness compared to the other methods. To elaborate, the Krylov subspace projection allows for an a posteriori error control mechanism based on error estimates that will be described below. On the other hand, polynomial expansion methods such as Taylor series or Chebyshev polynomial expansion offer no such a posteriori mechanism. While one can a priori estimate, based on the spectral radius ofÃ n , the number of terms required in the polynomial expansion for a desired accuracy, they remain highly over-estimated. Furthermore, in the case of split-operator, the efficacy of it rests on operating back and forth between Fourier and real space, so as to diagonalize the kinetic and the potential part of the Kohn-Sham Hamiltonian in succession. Thus, it involves the use of fast Fourier transforms (FFTs) which pose serious challenges to the parallel scalability of the code. The Krylov subspace projection, on the other hand involves operations only in real space and affords good parallel scalability (as will be shown in Sec. 5.4).
We now discuss the details of the Krylov subspace projection approach for the second-order Magnus operator. To begin with, a k-dimensional Krylov subspace for the matrixÃ n and the vectorψ is given by The Lanczos iteration provides a recipe for generating an orthonormal set of vectors Q k = {q 1 , q 2 , . . . , q k }, with q 1 =ψ/ ψ , that spans the same space as K k (Ã n ,ψ), and has the following property [79]: where T k is a tridiagonal matrix (upper Hessenberg matrix ifÃ n is non-Hermitian) and represents the projection ofÃ n into Q k (i.e., T k = Q T kÃn Q k ), e i is the i-th unit vector in C k , and β k+1,k is the (k + 1, k) entry of T k+1 = Q T k+1Ãn Q k+1 . By construction, we haveψ = ψ Q k e 1 . Furthermore, the operator Q k Q T k denotes the orthogonal projector onto the subspace K k (Ã n , ψ). Thus, one can obtain the best approximation to eÃ n ψ in the subspace K k (Ã n , ψ), denoted by y k , through a Galerkin projection, and is given by However, as evident, the above expression requires the evaluation of the action of eÃ n on Q k which remains computationally impractical beyond small system sizes, owing to the large size of the matrixÃ n . To this end, one can use a different approximation to eÃ n ψ, denoted by z k ∈ K k (Ã n , ψ), given by As is evident from the above form, the problem is now reduced to the evaluation of exp(T k ), wherein T k is a small matrix of size k × k, and hence, exp(T k ) can be evaluated inexpensively either through Taylor series expansion or exact eigenvalue decomposition of T k . The error, k , incurred in the above approximation is given by [80] Thus, the above relation provides a robust and inexpensive scheme to adaptively determine the dimension of the Krylov subspace by checking if k is below a set tolerance. An economic choice for the tolerance for k is determined from atomic RT-TDDFT calculations, such that it achieves < 10 meV accuracy in the excitation energies. For a multi-atom system, we employ the lowest such tolerance obtained for each of the constituent atomic species. Finally, we comment upon the numerical details of the second-order Magnus propagator with midpoint integration rule. As discussed in Sec. 3.3, the use of second-order Magnus propagator with midpoint integration rule, i.e., eÃ n , requires the knowledge ofH at a future time instant i.e.,H[t n−1 + ∆t/2], which is a priori unknown. A fully consistent approach involves, for a givenψ n−1 , the following steps: (i) approximatē arising out of the Lanczos procedure (Eq. 111) at each iterate of the self-consistent iteration. An efficient and sufficiently accurate approach is to use a predictor-corrector method to, first, predictH[t n−1 +∆t/4] through an extrapolation (linear or higher-order) from previous instants ofH, use it to propagateψ n−1 toψ n−1/2 , which is then used to evaluateH[t n−1 + ∆t/2]. We refer to [81] for the details of the predictor-corrector scheme. We remark that this predictor-corrector scheme is accurate to O(∆t 2 ), and hence, does not affect the estimates obtained in our full-discrete error analysis.

Results
In this section, we discuss the accuracy, rate of convergence, computational efficiency and the parallel scalability of higher-order finite-element discretization in conjunction with second-order Magnus propagator, for both pseudopotential and all-electron RT-TDDFT calculations. Based on the system, we use hexahedral spectral finite-elements of polynomial order 1 to 5, denoted as HEX8, HEX27, HEX64SPEC, HEX125SPEC, and HEX216SPEC. For the pseudopotential calculations, we provide comparison, in terms of accuracy and performance, of the higher-order finite-elements against the finite-difference method. The finite-difference based calculations are performed using the Octopus [30] software package. In all our finite-difference based calculations, we have used a stencil of order 4 in each direction, which is the default stencil order in Octopus. All the pseudopotential calculations are done using the norm-conserving Troullier-Martins pseudopotentials [82]. For all calculations, the ground-state Kohn-Sham orbitals are used as the initial states. We use the Chebyshev filter acceleration technique (see [45]) to efficiently compute the ground-state, for all the calculations done using finite-elements. All our scalability as well as benchmark studies demonstrating the computational efficiency are conducted on a parallel computing cluster with the following configuration: Intel Xeon Platinum 8160 (Skylake) CPU nodes with 48 processors (cores) per node, 192 GB memory per node, and Infiniband networking between all nodes for fast MPI communications.

Rates of Convergence
In this section, we study the rates of convergence of the dipole moment with respect to both finite-element mesh-size, h, as well as time-step, ∆t. We use methane and lithium hydride molecules as our benchmark systems for studying the rates of convergence, for pseudopotential and all-electron calculations, respectively.
We note that in order to study the convergence with respect to mesh-size, the dominant error must arise from spatial-discretization. To this end, we contain other sources of error, namely, time-discretization error and Krylov subspace projection error, by choosing a very small time-step of ∆t = 10 −4 , and using a small tolerance of 10 −12 for the Krylov subspace error (see Eq. 114). In effect, we mimic a semi-discrete (discrete in space but continuous in time) error analysis. Furthermore, to determine the coarsening rate of the mesh, as per Eq. 108, we use the a priori knowledge of the asymptotics of the ground-state atomic orbitals and electrostatic potentials. We remark that although we propose for using the far-field asymptotics of the ground-state electronic fields to determine the mesh coarsening rate, it, nevertheless, forms a reliable and cost effective way for discretizing the mesh away from the nuclei as opposed to any ad-hoc coarsening or using uniform discretization.
We now briefly touch upon on how to inexpensively obtain the single-atom ground-state electronic fields in the all-electron case. We note that for the single-atom case, under the assumption of equal fractional occupancy for degenerate orbitals, the electron density, and hence, the Kohn-Sham potential is spherically symmetric. Thus, one can rewrite the Kohn-Sham eigenvalue problem in spherical coordinates and invoke separation of variables to express a Kohn-Sham orbital as In the above equation, R nl (r) represents the radial part of the Kohn-Sham orbital corresponding to the principal quantum number n and azimuthal quantum number l, Y lm (θ, φ) denotes the real form of spherical harmonics for the pair of azimuthal quantum number l and magnetic quantum number m, and θ and φ represent the polar and azimuthal angles, respectively. The radial part, R nl (r), can be solved using the following 1D radial eigenvalue equation where nl is the eigenvalue corresponding to R nl (r). Similarly, the spherically symmetric electrostatic potentials, V H (r) and V ae N (r), can be obtained from the following 1D radial equations where Z denotes the atomic number of the atom under consideration and δ(r) is the Dirac-delta distribution. We note that all the above radial equations can be pre-computed inexpensively using 1D finite-elements, and stored for various atomic species for future use in obtaining optimal mesh coarsening rates. Subsequently, for a multi-atom case, we use the superposition of atomic orbitals and electrostatic potentials, in Eq. 108, to determine the coarsening rate. The resultant mesh forms the coarsest mesh for our convergence studies with respect to mesh-size. We then refine the coarsest mesh near the nuclei to obtain a sequence of increasingly refined meshes. We employ finite-elements of three different orders-HEX8, HEX27, and HEX125SPEC-in all our convergence studies. We remark that, although in the above procedure the single electronic fields are obtained for all-electron case, the same single-atom fields can be used for obtaining the coarsening rates for pseudopotential calculations as well. This is owing to the fact that the pseudo-wavefunctions and the pseudo-electrostatic potentials have the same far-field decay rates as its all-electron counterparts. Alternatively, one can directly use the single-atom pseudopotentials and pseudo-wavefunctions, obtained from a pseudopotential generation program, like FHI98PP [83].
We note that in order to perform the convergence study with respect to mesh-size, and compare to the semi-discrete error estimate obtained in Eq. 105, we require the knowledge of the continuous value of the dipole moment, µ x (t). To this end, we use the discrete dipole moment µ h x (t), obtained from a sequence of increasingly refined HEX125SPEC finite-element meshes to obtain a least-square fit of the form to determine µ x,0 (t), C, and q. In the above equation h min represents the minimum element size, h e , present in the mesh. The obtained µ x,0 (t) represents the extrapolated continuum limit (continuous in space) for the dipole moment computed using HEX125SPEC element, and is used as the reference value to compute |µx,0(t)−µ h for HEX8 and HEX27 finite-elements.
Next, we consider the convergence with respect to temporal discretization, i.e., ∆t. To this end, we use a sufficiently refined HEX125SPEC finite-element mesh and use increasingly refined ∆t to obtain a least-square fit of the form to determine µ h x (t n ), C, and q. The value of µ h x (t n ) obtained from the above equation represents the extrapolated continuum limit (continuous in time) for the dipole moment at t n .

All-electron calculations: Lithium Hydride
In this example, we conduct all-electron RT-TDDFT study on a lithium hydride molecule (LiH) with Li-H bond-length of 3.014 a.u. . A large cubical domain of length of 50 a.u. is chosen to ensure that the electron density decays to zero on the domain boundary, thereby, allowing us to impose Dirichlet boundary condition on the time-dependent Kohn-Sham orbitals and the Hartree potential. We use the aforementioned mesh-coarsening strategy based on the knowledge of ground-state atomic electronic fields to adaptively coarsen the mesh away from the nuclei. We remark that this adaptive coarsening of the mesh away from the nuclei allows us to use large computational domain sizes without significantly increasing the number of elements. This, in turn, allows us to circumvent the need for artificial absorbing boundary conditions, which are otherwise essential to tackle wave reflection effects that are observed while dealing with small computational domains. For all the meshes under consideration, we first obtain the ground-state and employ a weak delta-kick to excite the system. To elaborate, we use an electric field of the form E 0 (t) = κδ(t)x, with κ = 10 −3 a.u., where δ(t) is the Dirac-delta distribution andx is the unit vector along x-axis. This amounts to perturbing the groundstate Kohn-Sham orbitals, ψ GS α , by a factor e −iκx . Thus, our initial-states are defined as ψ α (0) = e −iκx ψ GS α . Figure 1 depicts the rates of convergence for the dipole moment at t = 1.0 a.u. for different orders of finite-elements. For all the three types of finite-elements under consideration, we observe close to optimal rates of convergence, O(h p ), where p is the degree of the finite-element interpolating polynomial. As is evident from Figure 1, much higher accuracies are obtained with HEX125SPEC when compared to HEX8 and HEX27 of the same mesh size. In particular, 200, 000 HEX8 elements are required to achieve relative errors of 10 −2 , whereas we achieve relative error of 10 −3 with just 3, 000 HEX125SPEC elements. Next, we study the rate of convergence of the dipole moment with respect to temporal discretization. To this end, we used a sufficiently refined HEX125SPEC mesh which affords 10 −4 relative error with respect to spatial discretization. We then propagate the initial-states using second-order Magnus propagator with different ∆t. Figure 2 depicts the rate of convergence of the dipole moment with respect to ∆t at t n = 1.0 a.u.. We obtain the rate of convergence, q = 1.96 (defined in Eq. 119), which agrees remarkably well with the theoretical rate of convergence for second-order Magnus propagator (see Proposition 3.2).

Pseudopotential calculations: Methane (CH 4 )
We now turn to examining rates of convergence for the pseudopotential case. We use a methane molecule with C-H bond-length of 2.07846 a.u. and a C-H-C tetrahedral angle of 109.4712 • as our benchmark system. As was the case with lithium hydride, we use ground-state single-atom electronic fields to coarsen the mesh away from the nuclei. We, once again, make use of a large cubical domain of length 50 a.u. to mimic simulations in R 3 . For all the meshes, we first, obtain the ground-state and then excite the system using a Gaussian electric field of the form E 0 (t) = κe (t−t0) 2 /w 2x , with κ = 2 × 10 −5 a.u., t 0 = 3.0 a.u., and w = 0.2 a.u. . Figure 3 illustrates the rates of convergence of the dipole moment at t = 5.0 a.u. for different orders of finite-elements. As in the case of lithium hydride, we obtain close to optimal rates of convergence, and observe significantly higher accuracies for HEX125SPEC over HEX8 and HEX27. Next, we study the rate of convergence afforded by the second-order Magnus propagator with respect to the time-step. To this end, we use a sufficiently refined HEX125SPEC mesh and employ the second-order Magnus propagator with different ∆t to propagate the ground-state Kohn-Sham orbitals under the influence of the Gaussian electric field. Figure 4 shows the rate of convergence of the dipole moment with respect to ∆t at t n = 5.0 a.u. . As was the case with lithium hydride, we obtain a convergence rate that is remarkably close to optimal rate of convergence, with q = 1.98 (defined in Eq. 119).

Computational Cost
In this section, we investigate the relative computational efficiency afforded by higher-order finite-elements over linear finite-element. We consider the previous two systems, lithium hydride and methane, for allelectron and pseudopotential calculations, respectively. We use the same strategy of mesh coarsening as detailed in Sec. 5.1. Since the objective of this study is to compare the relative performance of various orders of finite-elements, we eliminate any time-discretization effect by setting ∆t = 10 −4 . Furthermore, we use a tolerance of 10 −12 for the adaptive Lanczos (see Eq. 114) in order to eliminate any Krylov subspace projection error influencing the spatial discretization error. We repeat the previous numerical studies by exciting the lithium hydride molecule with a delta-kick (see Sec. 5.1.1), and the methane molecule with a Gaussian electric field (see Sec. 5.1.2). Figures 5 and 6 show the relative error in the dipole moment against the normalized computational time, for three different orders of finite-elements. The normalization of the computational time is done with respect to the longest time among the various meshes under consideration. As is evident, the relative computational efficiency afforded by higher-order finite-elements improves as the desired accuracy is increased. In particular, for a relative accuracy of 10 −3 , HEX125SPEC outperforms HEX8 and HEX27 by factor 150 − 200 and 10 − 18, respectively. This underscores the efficacy of higherorder finite-elements for RT-TDDFT calculations, an aspect which had, heretofore, remained unexplored for RT-TDDFT.

Other materials systems
In this section, we investigate the accuracy and computational efficiency afforded by higher-order finiteelements for other materials systems, in both pseudopotential and all-electron RT-TDDFT calculations. We also provide comparison, wherever possible, against pseudopotential calculations based on a finite-difference discretization. We use Al 2 , Al 13 , and Mg 2 as the benchmark metallic systems for pseudopotential calculations. Furthermore, we use Buckministerfullerene (C 60 ) as our benchmark insulating system for pseudopotential calculations. For the all-electron case, we use methane and benzene as our benchmark systems. Additionally, for the all-electron calculations we provide comparison, in the absorption spectrum, with their pseudopotential counterparts. We remark that for all the above systems under consideration, except Mg 2 , we use weak electric fields to excite them. For Mg 2 , we use a strong laser pulse to study the efficacy of higher-order finite-elements for nonlinear response.

Pseudopotential calculations: Al 2
We consider an aluminum dimer (Al 2 ) of bond-length 4.74 a.u. . In order to generate a suitable mesh, we use HEX64SPEC finite-elements with local refinement around the nuclei and coarsening away from the nuclei. The coarsening, as in the case of previous examples, was done by using the superposition of groundstate atomic electronic fields in Eq. 108. The local refinement around the atom is done such that it achieves an accuracy of 10 meV in the ground-state energy per atom. We use a cubical domain of length 60 a.u. to ensure that the wavefunctions decay to zero, and thereby, avoid any reflection effects. We excite the system from its ground-state by using a Gaussian electric of the form E 0 (t) = κe (t−t0) 2 /w 2x , with κ = 2 × 10 −5 a.u., t 0 = 3.0 a.u., and w = 0.2 a.u. . We use the second-order Magnus propagator with a time-step of ∆t = 0.05 a.u. = 0.001209 fs to propagate the Kohn-Sham orbitals and allow the system to evolve for 10 fs. We use a tolerance of 10 −8 for the Krylov subspace projection error (see Eq. 114). Next, we take the Fourier transform of the dipole moment to obtain the dynamic polarizability, α a,b (ω), where a is the index of the electric field's polarization direction and b is the index of the measurement direction of the dipole. Subsequently, we obtain the absorption cross-section tensor σ(ω) = 4πω c Im[α(ω)], where c is the speed of light. Finally, the dipole strength function (absorption spectrum), S(ω), is given by S(ω) = 1 3 Tr[σ(ω)]. The peaks in the absorption spectrum correspond to the excitation energies. In this study, we also assess the performance of higher-order finite-elements by comparing against finite-difference based method. To this end, we use the Octopus [30] software package to perform the finite-difference calculation. In order to obtain a suitable grid in Octopus, we first obtain the optimal grid-spacing and domain size that achieves an accuracy of 10 meV in the ground-state energy per atom, commensurate with the accuracy targeted in the finite-element discretization. Typically, one needs a larger domain for RT-TDDFT calculations than ground-state calculations, so as to avoid reflection at the domain boundaries. However, unlike the case in finite-elements, Octopus uses a uniform mesh, which makes the use of large domain sizes prohibitively expensive. A typical workaround is to use a smaller domain with absorbing boundary conditions. In order to better assess the effects of domain size for the finite-difference mesh, we use two cubical domains of sizes 38 a.u and 46 a.u, both with a grid-spacing of 0.2 a.u. . Furthermore, to understand the effect of absorbing boundary conditions, we perform an additional finite-difference calculation on the 38 a.u. mesh with a negative imaginary potential (NIP) near the boundaries. In particular, we use a potential of the following form  Fig. 7 compares absorption spectrum obtained from finite-elements against finite-difference. We have a used a Gaussian window of the form g(t) = e −αt 2 , with α = 0.005 a.u., in the Fourier transform of the dipole moment to artificially broaden the peaks, which is a standard practice in the computation of the absorption cross-section. As is evident from the figure, we get good agreement with the finite-difference based results for the domain size of 46 a.u. . The finite-difference calculation with domain size 38 a.u., with and without the absorbing boundary condition, provides qualitatively different results with two peaks around 5 eV. We attribute this discrepancy to possible reflection effects from the boundary, as a domain size of 38 a.u. may not be sufficient to avoid finite-domain size effects. Furthermore, comparing FD-38 and FD-38-ABS curves, we observe that the use of NIP based absorbing boundary condition, on its own, hardly improves the answer. This suggests that, for the system under consideration, one cannot rely, solely, on absorbing boundary conditions to avoid reflection effects, and hence, must use a larger domain. Table 1 compares the first two excitation peaks, the degrees of freedom and the total computational time for the finite-element and the finite-difference (46 a.u. domain size) based calculations. As is evident from the table, both finite-element and finite-difference based results agree to within 10 meV in the excitation energies. Furthermore, in terms of computational efficiency, we observe a ∼ 65-fold speedup for finite-elements over finite-difference. We remark that this superior efficiency for the finite-elements is largely attributed to fewer degrees of freedom that one can afford in finite-elements due to adaptive coarsening of the mesh, as opposed to a uniform mesh in finite-difference. We underline this by noting that while finite-difference requires over 12 million degrees of freedom, the finite-elements require only 31,411 degrees of freedom to attain similar accuracies.    13 We now consider a 13 atom aluminum cluster with an icosahedral symmetry. We use the same characteristic finite-element mesh as that of Al 2 but with a cubical domain of length 70 a.u., to avoid reflection effects. We repeat the same study as in the case of Al 2 . To elaborate, we use the same electric field, time-step, duration of propagation, and tolerance for the Krylov subspace projection as in the case of Al 2 .
We, once again, provide a comparative study against finite-difference based calculation by using a uniform cubical mesh of size 56 a.u. and grid-spacing 0.2 a.u. . Fig. 8 compares absorption spectrum obtained from finite-elements against finite-difference. We have used the same Gaussian window as in the case of Al 2 . As is evident from the figure, the peaks for both finite-element and finite-difference are in good agreement. Table 2 compares the first two excitation peaks, degrees of freedom, and the total computational time for the finite-element and the finite-difference based calculations. Both the methods agree to within 10 meV in the first two excitation energies. In terms of computational efficiency, the finite-elements attain a ∼ 8-fold savings in the computational time against finite-difference, once again, attributed to the fewer degrees of freedom in finite-elements owing to adaptive coarsening of the mesh. In particular, the finite-elements afford ∼ 30-fold fewer degrees of freedom as compared to finite-difference.

Pseudopotential calculations: Buckministerfullerene
In this example, we consider the Buckministerfullerene molecule, as shown in Figure 9. It comprises of 60 carbon atoms (240 electrons) packed into the shape of a buckyball. In order to generate a suitable finite-element mesh, we use HEX64SPEC finite-elements with the local refinement around the nuclei and coarsening away from the nuclei, governed by Eq. 108. We use a cubical domain of length 50 a.u. to eliminate any reflection effects from the boundaries. The ground-state is excited with a Gaussian electric field of the same form as used for Al 2 . We use a time-step of 0.05 a.u. and propagate the system for 10 fs. As with previous cases, we also assess the performance of higher-order finite-elements by comparing against finite-difference based method, as implemented in the Octopus package. We assess the effects of domain size for the finite-difference mesh, by using two cubical domains of sizes 30 and 36 a.u, both with a grid-spacing of 0.15 a.u. . Furthermore, we study the effect of absorbing boundary conditions by performing an additional finite-difference calculation on the 30 a.u. mesh with a negative imaginary potential (NIP) near the boundaries. We use an NIP of the same form as used in Al 2 (see Sec. 5.3.1), albeit with L = 14.0 a.u. . We denote these three finite-difference calculations, namely, with domain size 36 a.u., with domain size 30 a.u., and with domain size 30 a.u. along with NIP absorbing boundary condition as FD-36, FD-30, and FD-30-ABS, respectively. Figure 10 shows the absorption spectrum obtained from finite-element and the three different finite-difference based calculations. We have used a same Gaussian window of the form g(t) = e −αt 2 , with α = 0.01 a.u., to artificially broaden the peaks. As is evident from the figure, there is good agreement between the finite-element and FD-36 for all the excitation peaks. On the other hand, while FD-30 and FD-30-ABS have good agreement with finite-elements for the first two peaks, they differ for the rest, possibly because of reflection effects. Furthermore, comparing FD-30 and FD-30-ABS, we remark that the use of NIP based absorbing boundary condition is not improving the absorption spectrum. This, once again, indicates that one cannot always dispense using a larger domain by, solely, using absorbing boundary conditions. Table 3 compares the first two excitation peaks, degrees of freedom, and the computational time for finite-elements against that of FD-36. Both finite-element and FD-36 based results match within 30 meV in the first two peaks. In terms of computational efficiency, finite-elements attain a ∼ 3-fold speedup over FD-36. This higher efficiency of the finite-elements, is once again, attributed to a ∼ 9-fold fewer degrees of freedom required by the finite-elements against that of finite-difference.  In this example, we study the higher harmonic generation in a magnesium dimer with bond-length of 4.74 a.u. . Unlike the previous examples, we use a strong laser pulse to excite the system from its ground-state. In particular, we use an electric field of the form where E 0 , ω l , and T denote the amplitude, frequency, and duration of the laser pulse, respectively. In our example, we use E 0 = 0.01, ω l = 0.03, and T = 5 2π ω l (all in atomic units). We use HEX125SPEC finite-elements with the local refinement around the nuclei and coarsening away from the nuclei, governed by Eq. 108. Furthermore, we use a cubical domain of length 100 a.u. to eliminate any reflection effects from the boundaries. For time-propagation, we use the second-order Magnus propagator with ∆t = 0.025 and propagate it for 25.33 fs (i.e., T fs). We use a tolerance of 10 −8 for the Krylov subspace projection error (see Eq. 114). Next, we obtain the dipole power spectrum, P (ω), of the system by taking the imaginary part of the Fourier transform of the acceleration of the dipole moment, µ(t). To elaborate, Theoretically, for a system with spatial inversion symmetry, only odd multiples of the frequency of the exciting laser pulse must be emitted. We verify this in Figure 11 wherein the peaks in the power spectrum coincide with odd harmonics. Furthermore, we observe that the decay of the intensity of the peaks flattens beyond the 13-th harmonic, which corroborates well with the plateau phenomenon, typically observed in experiments [84]. We emphasize that despite the large domain size used in this calculation, we require only ∼ 60, 000 basis functions. This underlines the efficacy of higher-order finite-elements for even nonlinear regime in RT-TDDFT.

All-electron calculations: Methane (CH 4 )
We now examine the competence of higher-order finite-elements for all-electron RT-TDDFT calculations by providing a comparative study with its pseudopotential counterpart. In this example, we consider a methane molecule with the same geometry as described in Sec. 5.1.2. We use HEX64SPEC and HEX125SPEC elements for the pseudopotential and all-electron case, respectively. For both all-electron and pseudopotential cases, we use the same mesh coarsening strategy as used in all previous examples. For both the meshes, we use a large cubical domain of length 40 a.u., so as to eliminate reflection from the boundaries. Both the systems are excited from their respective ground-states using the same Gaussian electric field as used for Al 2 (see Sec. 5.3.1) and propagate them for 10 fs. We use time-steps of 0.05 and 0.025 (both in a.u.) for the pseudopotential and all-electron calculations, respectively. Furthermore, a tolerance of 10 −8 for the Krylov subspace projection error (see Eq. 114) is used for both the calculations. The absorption spectra for both the calculations are shown in Figure 12. We used the same Gaussian window as in the case of  5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 Power Spectrum Harmonic Order Buckministerfullerene (see Sec. 5.3.3), to artificially broaden the peaks. As evident from the figure, we obtain remarkable agreement between the all-electron and pseudopotential results, i.e., the two curves are almost identical. Table 4 we list the first two excitation peaks, degrees of freedom, and total computational time for both the calculations. The first two excitation peaks agree to within 10 meV. We remark that the all-electron calculation requires ∼ 100x more computational time as compared to the pseudopotential case. This large computational expense for the all-electron calculation stems primarily from the need of a highly refined mesh near the nuclei, so as to accurately capture the sharp variations in the electronic fields near the nuclei. This refinement has two major consequences: a) an increase in the degrees of freedom; and b) widening of the spectral radius of the Kohn-Sham Hamiltonian, which, in turn, warrants a smaller time-step (see Eq. 109) as well as a larger Krylov subspace to achieve the prescribed accuracy (see Eq. 114). In particular, for the case of methane, we required ∼ 4x degrees of freedom and ∼ 10x the size of the Krylov subspace as compared to that of the pseudopotential case. We emphasize that while finite-elements are expensive for the all-electron calculation, they provide the desired accuracy and offer systematic convergence (see Sec. 5.1.1). Moreover, one can mitigate the need of a refined mesh for the all-electron calculation by using an enriched finite-element basis, wherein the standard (classical) finite-element basis are augmented with numerical atom-centered basis. This idea has successfully attained 100−300x speedup over the standard (classical) finite-elements for ground-state DFT calculations [52], and can be extended to RT-TDDFT to further the capabilities of finite-elements.  In this example, we perform similar comparative studies between all-electron and pseudopotential calculations for benzene molecule. As with the methane molecule, we use HEX64SPEC and HEX125SPEC finite-elements for the pseudopotential and all-electron calculation, respectively. Furthermore, we use the same characteristic mesh features (i.e., refinement near the nuclei, coarsening rate, simulation domain), in both the meshes, as their counterparts in the methane calculation (see Sec. 5.3.5). Rest of the simulation details, namely, time-steps, duration of propagation, and the tolerance for the Krylov subspace projection error, remain the same as that used for methane. Figure 13 compares the absorption spectra from the allelectron and pseudopotential calculations. While there is qualitative agreement between both calculations, quantitatively the predictions from all-electron and pseudopotential calculations differ. In particular, the first two excitation peaks (see Table 5) differ up to ∼ 0.2 eV. This suggests that one ought to carefully test for the transferability of the pseudopotential approximation used, to provide reliable quantitative predictions from RT-TDDFT calculations. We take note that a more careful comparison of pseudopotential and all-electron calculations warrants a scan through a range of pseudopotential approximation. Nevertheless, the objective of this exercise is to highlight the fact that finite-elements, by treating both pseudopotential and all-electron calculations on an equal footing, allows for a robust tool for such transferability studies. Table 5: Comparison of all-electron (AE) and pseudopotential (PSP) calculations for benzene: First and second excitation energies (E 1 , E 2 , respectively, in eV), degrees of freedom (DoF), and total computation CPU time (in CPU hours).

Scalability
We demonstrate the parallel scalability (strong scaling) of the proposed enriched finite element basis in Figure 14. We choose the Buckministerfullerene molecule containing ∼ 3.5 million degrees of freedom AE PSP Figure 13: Absorption spectra of benzene (number of basis functions) as our fixed benchmark system and report the relative speedup with respect to the wall time on 24 processors. The use of any number of processors below 24 was unfeasible owing to the memory requirement posed by the system. As is evident from the figure, the scaling is in good agreement with the ideal linear scaling behavior up to 384 processors, at which we observe a parallel efficiency of 86.2%. However, we observe a deviation from linear scaling behavior at 768 processors with a parallel efficiency of 74.2%. This is attributed to the fact that, at 768 processors, the number of degrees of freedom possessed by each processor falls below 5000, which is low to achieve good parallel scalability.

Summary
In summary, in this work, we have investigated the accuracy, computational efficiency and scalability of higher-order finite-elements for the RT-TDDFT problem, for both pseudopotential and all-electron calculations. We begin by decomposing the discretization error in the time-dependent Kohn-Sham orbitals into two parts, one arising due to spatial discretization and the other due to temporal discretization. We, then, developed semi-discrete error estimates (discrete in space, continuous in time) in the Kohn-Sham orbitals in terms of the electronic fields, and the order and size of the finite-elements used. Subsequently, we developed full-discrete error-estimates in the Kohn-Sham orbitals, in the context of second-order Magnus propagator with a midpoint integration rule, in terms of the time-step size and the spatial derivative of the Kohn-Sham orbitals. We note that finite-elements form a non-orthogonal basis, and hence, require efficient strategies to numerically evaluate the inverse of the overlap matrix in order to construct the Magnus propagator (or any exponential time-integrator). To this end, we employed spectral finite-elements in conjunction with Gauss-Legendre-Lobatto quadrature rules that results in a diagonal overlap matrix, thereby making the evaluation of its inverse trivial. We note that the second-order Magnus propagator with a midpoint integration rule warrants the knowledge of the time-dependent Kohn-Sham Hamiltonian at half time-step in future, which is unknown a priori. To this end, we used a second-order predictor-corrector scheme to accurately evaluate the Hamiltonian at a future time instance. Furthermore, we employed an adaptive Lanczos subspace projection to evaluate the action of the Magnus propagator, defined as exponential of a matrix, on the Kohn-Sham orbitals.
In our numerical studies, we used the semi-discrete error estimates along with the a priori knowledge of the asymptotics of the electronic fields, using superposition of ground-state atomic solutions, to adaptively coarsen the mesh away from the nuclei and obtain close to optimal finite-element meshes, for both pseudopotential and all-electron systems. We used methane and lithium hydride molecule as our benchmark systems for pseudopotential and all-electron calculations, respectively. To begin with, we investigated the rates of convergence of the dipole moment with respect to mesh size for linear, quadratic and fourth-order spectral finite-elements. In both the benchmark systems, we observed close to optimal rates of convergence for all the three different orders of finite-elements employed. Furthermore, for a fixed order of finite-element, we investigated the rate of convergence of the dipole moment with respect to time-step as afforded by the second-order Magnus propagator. We attained close to optimal rates of convergence with respect to time discretization for both pseudopotential and all-electron calculations. Next, we assessed the computational efficiency afforded by higher-order finite-elements over linear finite-elements, using, once again, methane and lithium hydride as our benchmark systems for pseudopotential and all-electron calculations, respectively. To this end, we used the developed mesh coarsening scheme and studied the computational time required to solve the benchmark systems. We find a staggering 100−fold and 10−fold speedup afforded by the fourth-order finite-element over linear and quadratic finite-elements, respectively. Furthermore, we assessed the accuracy and efficiency afforded by higher-order finite-elements for other materials systems. We used Al 2 , Al 13 , and Buckministerfullerene (C 60 ) as the benchmark systems for pseudopotential calculations and compared the resulting absorption spectrum to that obtained from finite-difference based calculations done using the Octopus package. We observed good agreement in the excitation energies obtained from both the methods. In terms of computational efficiency, we obtained 40x, 8x, 3x speedup over finite-difference for Al 2 , Al 13 , and Buckministerfullerene, respectively. Thus, we underscore that higher-order finite-elements can compete with finite-difference, especially for non-periodic systems. Furthermore, we remark that, with the adaptive mesh coarsening afforded by finite-elements, one can use large simulation domains without significant increase in the number of basis functions, and hence, circumvent the use of any artificial absorbing boundary conditions. We also demonstrated the efficacy of higher-order finite-elements, especially its efficient handling of large domains, for nonlinear response by studying the higher harmonic generation in Mg 2 under a strong electric field. Next, we conducted similar studies for the all-electron case, using methane and benzene as our benchmark systems. For both the systems, we obtained good agreement with corresponding pseudopotential based results, thereby, establishing the capabilities of finite-elements for all-electron calculations. Additionally, the all-electron vis-á-vis pseudopotential results highlight the versatility of the finite-elements in handling both pseudopotential and all-electron calculations on an equal footing. Lastly, in terms of parallel scalability, we obtained good parallel efficiency with almost linear scaling up to 384 processors for a benchmark system comprising of the Buckministerfullerene molecule containing ∼ 3.5 million basis functions, and good scaling up to 768 processors.
Thus, higher-order finite-elements offer a computationally efficient, systematically improvable, and scalable basis for RT-TDDFT calculations, applicable to both pseudopotential and all-electron cases. We note that, while for all-electron calculations one needs a large number of basis functions in order to capture the sharp electronic fields near the nuclei, this shortcoming can be alleviated by augmenting the finite-element basis with numerical atom-centered basis. This idea has been successfully used for ground-state DFT [52], and its extension to RT-TDDFT is currently being investigated. Furthermore, the key ideas behind the a priori mesh adaptation developed in this work can be used to formulate an a posteriori mesh refinement scheme, thereby, allowing for more robust error control. Lastly, the prospect of using higher-order finiteelements to develop linear scaling RT-TDDFT algorithms is also of significance, and is a subject of our future investigation.