Stochastic Time-Dependent DFT with Optimally Tuned Range-Separated Hybrids: Application to Excitonic Effects in Large Phosphorene Sheets

TL;DR

This paper introduces a stochastic time-dependent DFT method with optimally tuned range-separated hybrids that efficiently calculates excitonic optical spectra in large phosphorene sheets, matching experimental results.

Contribution

It presents a novel stochastic approach to TDDFT with non-local exchange, enabling efficient and accurate excitonic spectra calculations for large 2D materials.

Findings

Accurately reproduces experimental optical spectra.

Scales quadratically with the number of electrons.

Demonstrates efficiency on large phosphorene sheets with up to 1958 electrons.

Abstract

We develop a stochastic approach to time-dependent DFT with optimally-tuned range-separated hybrids containing non-local exchange, for calculating optical spectra. The attractive electron-hole interaction, which leads to the formation of excitons, is included through a time-dependent linear-response technique with a non-local exchange interaction which is computed very efficiently through a stochastic scheme. The method is inexpensive and scales quadratically with the number of electrons, at almost the same (low) cost of time dependent Kohn-Sham (TDKS) with local functionals. Our results are in excellent agreement with experimental data and the efficiency of the approach is demonstrated on large finite phosphorene sheets containing up to 1958 valence electrons.

Figures

Equations

$$H[n,\rho ]\equiv h_0+v_H\left[n\left(\mathbf{r}\right)\right]+v_C^{\gamma }\left[n\left(\mathbf{r}\right)\right]+X^{\gamma }\left[\rho \left(\mathbf{r},{\mathbf{r}}^{\prime }\right)\right],$$

$$\frac{1}{r}=\frac{\mathrm{erfc}\left(\gamma r\right)}{r}+\frac{\mathrm{erf}\left(\gamma r\right)}{r},$$

$$X_{\mathrm{nl}}^{\gamma }(\mathbf{r},{\mathbf{r}}^{\prime })\equiv -\nu (\mathbf{r},{\mathbf{r}}^{\prime })\rho \left(\mathbf{r},{\mathbf{r}}^{\prime }\right).$$

$$⟨\mathbf{r}\left|X_{\mathrm{nl}}^{\gamma }\right|\psi ⟩=$$

$$-\sum _i{f}_i{\varphi }_i(\mathbf{r})\int \nu (\mathbf{r},{\mathbf{r}}^{\prime }){\varphi }_i\left({\mathbf{r}}^{\prime }\right)\psi \left({\mathbf{r}}^{\prime }\right)\mathrm{d}{\mathbf{r}}^{\prime }.$$

$$i|{\dot{\varphi }}_i\left(t\right)⟩=[H\left[n\left(t\right),\rho \left(t\right)\right]+\delta \nu (\mathbf{r},t)]|{\varphi }_i\left(t\right)⟩.$$

$$\eta \left(\mathbf{r}\right)=\sum _i\sqrt{f_i}{\varphi }_i\left(\mathbf{r}\right)⟨{\varphi }_i|\bar{\eta }⟩,$$

$$\rho ={\left\{|\eta ⟩⟨\eta |\right\}}_{\bar{\eta }}$$

$$\frac{\mathrm{erf}\left(\gamma \left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\right)}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}={\left\{\zeta \left(\mathbf{r}\right)\zeta \left({\mathbf{r}}^{\prime }\right)\right\}}_{\bar{\theta }}$$

$$\zeta \left(\mathbf{r}\right)=\frac{1}{V^{\frac{1}{2}}}\sum _{\mathbf{k}}{}^{{}^{\prime }}\sqrt{\nu \left(\gamma ,\mathbf{k}\right)}{e}^{i\left(\bar{\theta }\left(\mathbf{k}\right)+\mathbf{k}\cdot \mathbf{r}\right)}$$

$$⟨\mathbf{r}\left|X_{\mathrm{nl}}^{\gamma }\right|\varphi ⟩=$$

$$-{\left\{\xi (\mathbf{r})⟨\xi |\varphi ⟩\right\}}_{\bar{\eta }\theta }\simeq -\frac{1}{N_{\xi }}\sum _{j=1}^{N_{\xi }}{\xi }^j(\mathbf{r})⟨{\xi }^j|\varphi ⟩$$

$$\beta \left(\mathbf{r},t\right)=\sum _j{e}^{i{\alpha }_j(t)}\sqrt{f_j}{\varphi }_j\left(\mathbf{r},t\right),$$

$$⟨\mathbf{r}|X_{\mathrm{nl}}^{\gamma }(t)|\varphi ⟩$$

$$\times \int v(\mathbf{r}-{\mathbf{r}}^{\prime }){\beta }^{l\ast }({\mathbf{r}}^{\prime },t)\varphi ({\mathbf{r}}^{\prime })d{\mathbf{r}}^{\prime }.$$

$$i|{\dot{\varphi }}_i^{\Delta }\left(t\right)⟩$$

$$i|{\dot{\varphi }}_i^{\Delta =0}\left(t\right)⟩$$

$$\mu \left(t\right)=\frac{1}{\Delta }\int (\mathbf{r}\cdot \hat{\mathrm{e}})\left[n^{\Delta }\left(\mathbf{r},t\right)-n^{\Delta =0}\left(\mathbf{r},t\right)\right]\mathrm{d}\mathbf{r},$$

$${\varphi }_i(\mathbf{r},t=0^{+})=e^{-i(\mathbf{r}\cdot \hat{\mathrm{e}})\Delta }{\varphi }_i(\mathbf{r})$$

$$|{\varphi }_i(t+$$

$$e^{-i{X}_{\mathrm{nl}}^{\gamma }\frac{dt}{2}}$$

$$e^{-i{X}_{\mathrm{nl}}^{\gamma }\frac{dt}{2}}|\varphi ⟩\simeq N_{\varphi }\left(1-i{X}_{\mathrm{nl}}^{\gamma }\frac{dt}{2}\right)|\varphi ⟩$$

$$E_b=E_g-E_1,$$

$$\sigma \left(\omega \right)=$$

$$\frac{4\pi }{c}\omega \iint \delta \tilde{v}\left(\mathbf{r},\omega \right)\tilde{\chi }\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\delta \tilde{v}\left({\mathbf{r}}^{\prime },\omega \right)\mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime },$$

$$\chi \left(\mathbf{r},{\mathbf{r}}^{\prime },t-t^{\prime }\right)=\frac{\delta n\left(\mathbf{r},t\right)}{\delta v\left({\mathbf{r}}^{\prime },t^{\prime }\right)},$$

$$\sigma \left(\omega \right)=\frac{4\pi }{c}\omega \iint (\hat{\mathbf{e}}\cdot \mathbf{r})\tilde{\chi }\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\cdot (\hat{\mathbf{e}}\cdot {\mathbf{r}}^{\prime })\mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }.$$

$$\delta \nu ({\mathbf{r}}^{\prime },t)=(\hat{\mathbf{e}}\cdot {\mathbf{r}}^{\prime })\delta \left(t\right)\Delta ,$$

$$\mu (t)=\frac{1}{\Delta }\int (\hat{\mathbf{e}}\cdot \mathbf{r})\cdot \delta n\left(\mathbf{r},t\right)\mathrm{d}\mathbf{r}.$$

$$\sigma \left(\omega \right)=\frac{4\pi \omega }{c}{\int }_0^{\infty }\mu (t)e^{i\omega t}\mathrm{d}t.$$

Full text

I Introduction

The understanding of optical excitations in materials is essential for developing novel optical and electronic devices.Onida, Reining, and Rubio (2002); Marques and Gross (2004) However, it is very challenging to calculate optical properties in extended systems. For small molecules, highly correlated approaches are used, including configuration interaction or the equation of motion coupled cluster approach.Bartlett and Musiał (2007); Helgaker, Jorgensen, and Olsen (2014) Further, the optical response is typically described by the Bethe-Salpeter equation.Onida, Reining, and Rubio (2002); Martin, Reining, and Ceperley (2016) These high level techniques are predictive but scale steeply with the number of electrons so they can only be used for relatively small molecules and unit cells.

An alternative to costly many-body techniques is time-dependent density functional theory (TDDFT)Runge and Gross (1984) that describes excited state energies, geometries, and other properties of small molecules with a relatively moderate computational cost. In principle TDDFT is exact, but in practice approximations have to be introduced. The most common is the adiabatic time-dependent Kohn-Sham theory (TDKS). While TDKS has been applied successfully to a wide range of molecular systems,Onida, Reining, and Rubio (2002); Marques and Gross (2004) it suffers from many failures, particularly for extended systems, charge-transfer excited states,Dreuw, Weisman, and Head-Gordon (2003) multiple excitations,Maitra et al. (2004) and avoided crossings.Levine et al. (2006) The most notable problem of TDKS is the inability to capture low-lying excitonic states in bulk.Onida, Reining, and Rubio (2002); Martin, Reining, and Ceperley (2016)

It has been argued that a TDDFT formulation beyond the Kohn-Sham picture, namely, a TD-GKS (Generalized Kohn-Sham) approach Baer and Kronik (2018) which employs a non-local exchange interaction,Seidl et al. (1996); Baer and Neuhauser (2005); Kümmel and Kronik (2008) captures the necessary physics to describe excitation in extended systems Yang, Sottile, and Ullrich (2015); Refaely-Abramson et al. (2015); Brawand et al. (2016) and accurately predicts the formation of bound excitons. However, the inclusion of the non-local exchange in the TDDFT calculations makes them computationally demanding and out of reach for large nanoscale systems. Recently a family of stochastic orbital methods has been developed to describe ground and excited states with the goal of lowering the computational complexity at the cost of introducing a controllable statistical error. Baer, Neuhauser, and Rabani (2013); Neuhauser et al. (2014); Gao et al. (2015); Neuhauser et al. (2015); Rabani, Baer, and Neuhauser (2015); Vlček et al. (2017) Specifically relevant to the current work is a stochastic ground state range-separated hybrid DFT method Neuhauser et al. (2015) and a stochastic method for the Bethe-Salpeter equation (BSE). Rabani, Baer, and Neuhauser (2015)

Our first and main aim in this work is to overcome, using stochastic methods, the computational bottleneck in TDDFT with non-local exchange. Specifically, we develop a real-time generalized Kohn-Sham method based on a range-separated hybrid (RSH) with a long-range exact exchange operator. The approach has similarities to that proposed for stochastic BSE (where a damped exchange operator was used), but the starting point is different. Here, the starting point is long-range-corrected RSH-DFT, a method which is known to produce excellent charge-transfer states.Stein, Kronik, and Baer (2009) Therefore the present long-range exchange TDDFT starts only from a long-range DFT calculation (which is also implemented stochastically), so the resulting approach is self-contained in the DFT/TDDFT framework and does not resort to a separate computation of individual quasiparticle states; this is in contrast to a BSE work where the starting point is a prior calculation of the quasiparticle states which is achievable, e.g., through the stochastic $GW$ method.Neuhauser et al. (2014); Vlček et al. (2017, 2018)

The second aim of the paper is then to use the resulting fully ab-initio stochastic TDDFT method for describing optical excitations in extremely large systems with thousands of electrons.

Below, we first (Sec. II) review the basic theory and present our stochastic implementation. In Sec. III we verify our method by comparing with experiment for ${\rm PH_{3}}$ and the method is then applied to study 2D phosphorene sheets of increasing sizes. Conclusions follow in Sec. IV.

II Theory

It is well-known that absorption can be determined, in linear response, from the time-evolution of an induced dipole density (see Appendix I). This time evolution is in principle governed by the time-dependent Schrödinger equation, which is however intractable beyond few-electron systems. DFTHohenberg and Kohn (1964) is practical approach for recasting the many-electron system as a set of virtual particles interacting via a mean-field exchange-correlation (xc) potential.

The original formulation of the Kohn-Sham (KS) Kohn and Sham (1965) DFT scheme describes the xc interactions by a local potential which is in principle non-local in time. In practice it is further approximated, e.g., by (semi)local functionals of the electronic density at given time (i.e., the functional is adiabatic). As mentioned in the introduction, this formulation has some notable failures, e.g., missing excitonic effects. An alternative route, which we pursue here, is to employ a GKS schemeSeidl et al. (1996); Savin and Flad (1995); Baer and Kronik (2018) with non-local long-range asymptotic behavior as required for correct description of charge transfer and polarizability Ghosez, Gonze, and Godby (1997) and for electron-hole bound states.Baer and Neuhauser (2005); Yang, Sottile, and Ullrich (2015); Refaely-Abramson et al. (2015); Brawand et al. (2016)

We first review below the formulation of DFT and TDDFT with long-range non-local exchange, followed by the details of stochastic implementation.

II.1 DFT with long-range non-local exchange

The GKS Hamiltonian reads

[TABLE]

where $h_{0}$ contains the kinetic energy and the electron-nuclear attraction. The density-density repulsion is given by the Hartree potential $v_{H}$, and $v_{C}^{\gamma}$ (where $\gamma$ is defined shortly) is a (semi)local correlation density functional – we use here a local functional form.Perdew and Wang (1992) The non-local exchange interaction $X^{\gamma}$ is a functional of the density-matrix $\rho\left({\bf r},{\bf r^{\prime}}\right)$, where the density is of course $n({\bf r})=\rho({\bf r},{\bf r}).$ (Note that we use different symbols, $n$ and $\rho$, for the density and density matrix since later we calculate the two separately - one deterministically and the other stochastically.)

The form of $X^{\gamma}$ derives from a screened Coulomb interaction governed by a single parameter $\gamma$.Savin and Flad (1995); Leininger et al. (1997); Baer and Neuhauser (2005) Specifically, the Coulomb kernel is partitioned as:

[TABLE]

where $\gamma$ is the range separation parameter. The first term dominates at small distances ($r\to 0$) and its contribution to the exchange is approximated by a local density functional.Perdew and Wang (1992); Leininger et al. (1997); Livshits and Baer (2007) The second term in Eq. (2) is active at large distances and gives the non-local Fock exchange ($X_{{\rm nl}}^{\gamma}$). The matrix element of the exchange vector is a direct product of the density matrix and the non-local Coulomb interaction (the simple product form is important later in the stochastic formulation):

[TABLE]

Here, the long-range interaction is $\nu({\bf r},{\bf r^{\prime}})={\rm erf}\left(\gamma|{\bf r}-{\bf r^{\prime}}|\right)/|{\bf r}-{\bf r^{\prime}}|$, and the density matrix is determined from the eigenstates $\rho({\bf r},{\bf r^{\prime}})=\sum_{i}f_{i}\phi_{i}({\bf r})\phi_{i}({\bf r^{\prime}})$, where $i$ is a state and spin index and $f_{i}$ are the occupation factors, while the density is $n({\bf r})=\rho({\bf r},{\bf r})=\sum_{i}f_{i}|\phi_{i}({\bf r})|^{2}$. (In the following we do not denote spin explicitly.) Thus, deterministically, the nonlocal exchange term acts on a general function $\psi$ as

[TABLE]

In the first (DFT) stage, the occupied eigenstates $\phi_{i}({\bf r})$ of Eq. (1) are calculated self-consistently, ensuring $H\phi_{i}=\varepsilon_{i}\phi_{i},$ where the density and density matrix are functions of the eigenstates. The value of $\gamma$ is found by enforcing the IP theorem that requires that the HOMO energy equals the ionization energy. This optimal tuning leads to good IPs and fundamental band gaps in finite systems.Stein et al. (2010, 2012); Kronik et al. (2012)

The method’s success stems from the combination of (semi)local functionals that capture correlation effects well at short distances with the nonlocal $X_{{\rm nl}}^{\gamma}$ that guarantees for finite systems the asymptotically correct $1/r$ behavior of the exchange potential which is crucial for proper inclusion of the attractive electron-hole interaction.Ghosez, Gonze, and Godby (1997); Onida, Reining, and Rubio (2002)

II.2 TDDFT with long-range non-local exchange

It is well-known (Appendix I) that the absorption spectrum is obtained from a linear-response propagation of the density. Specifically, for polarized excitation along a unit vector $\hat{\boldsymbol{{\rm e}}}$, we apply a small perturbation $\delta\nu({\bf r},t)=({\bf r}\cdot\hat{\boldsymbol{{\rm e}}})\delta(t)\Delta$ where $\Delta$ is a small constant (typically between $10^{-3}$and $10^{-5}$ a.u.). Then, the system evolves under the time-dependent GKS equation (using $\hbar=1$)

[TABLE]

To simplify the notation we usually do not denote the dependence of the density and density matrix (and therefore of the time-dependent Hamiltonian) on the excitation strength $\Delta$.

The Hamiltonian $H\left(t\right)$ is time dependent as it explicitly depends on the propagated eigenstates $\phi_{i}\left({\bf r},t\right)$, the time-dependent charge density $n\left({\bf r},t\right)=\sum_{i}f_{i}|\phi\left({\bf r^{\prime}},t\right)|^{2}$, and the charge density matrix, $\rho\left({\bf r},{\bf r^{\prime}},t\right)=\sum_{i}f_{i}\phi\left({\bf r},t\right)\phi_{i}^{*}\left({\bf r^{\prime}},t\right)$. The dipole moment along the excitation direction is then calculated from the density, $\mu(t)=\int({\bf r}\cdot\hat{\boldsymbol{{\rm e}}})n({\bf r},t)d{\bf r},$ and the absorption spectrum is calculated by Fourier transforming the dipole moment $\mu(t)$ (Appendix I).

In principle, the exchange-correlation term in the Hamiltonian should account for memory effects, but since its form is unknown, we resort to the adiabatic approximation and construct the xc terms directly from $n\left({\bf r},t\right)$ and $\rho\left({\bf r},{\bf r^{\prime}},t\right)$. Thus, the difference from Kohn-Sham type adiabatic TDDFT is only in the exchange kernel.

The application of the non-local exchange as presented in Eq. (3) is computationally demanding, due to the integral over the density matrix. Practical calculations are therefore limited to systems with a low number of states.

II.3 Stochastic DFT with non-local exchange

Next we review our implementation (Neuhauser et al., 2015) of the DFT equations with a stochastic representation of the non-local exchange operator. This is followed by implementation of stochastic TDDFT in the next section.

The first step is the DFT ground-state calculation, where we use the stochastic-exchange approach of Ref. Neuhauser et al., 2015. This grid-based method is done by two key parts. The first is the representation of the density matrix as an average over stochastic correlation functions. Specifically, we construct stochastic states $\left\{\eta\right\}$, each of which is a linear combination of all the occupied eigenstates $\left\{\phi_{i}\right\}$ (cf. Refs.Baer and Neuhauser, 2012; Neuhauser, Rabani, and Baer, 2013; Baer, Neuhauser, and Rabani, 2013; Gao et al., 2015; Rabani, Baer, and Neuhauser, 2015)

[TABLE]

where $\bar{\eta}$ is a completely random real vector, e.g., $\bar{\eta}({\bf r})=\pm(dV)^{-1/2},$ and $dV$ is the grid volume-element. It is straightforward to show that as an operator, the density matrix becomes an average over the separable terms

[TABLE]

i.e., $\rho\left({\bf r},{\bf r^{\prime}}\right)=\left\{\eta\left({\bf r}\right)\eta\left({\bf r^{\prime}}\right)\right\}_{\bar{\eta}},$ where $\left\{\cdots\right\}_{\bar{\eta}}$ denotes a statistical average over all random states $\bar{\eta.}$ Since the average of $|\eta\rangle\langle\eta|$ yields the density matrix, we can view $\eta({\bf r})$ as a stochastic density amplitude.

In the ground-state DFT stage we supplement the stochastic representation of the density matrix by a similar stochastic decomposition of the long-range Coulomb interaction (Eq. (7)) using stochastic states $\zeta$:

[TABLE]

that are evaluated as

[TABLE]

where $\bar{\theta}\left({\bf k}\right)$ is a random phase and we impose $\bar{\theta}(-{\bf k})=-\bar{\theta}({\bf k})$ to ensure that $\zeta({\bf r})$ are real. Also, $V$ is the total volume. The prime in the summation indicates that the ${\bf k}=0$ term is excluded and is later added analytically. The $\gamma$-dependent long-range Coulomb interaction in momentum space is $\nu\left(\gamma,{\bf k}\right)=4\pi e^{-\frac{|{\bf k}|^{2}\gamma^{2}}{4}}/|{\bf k}|^{2}$. Note that the average in Eq. (8) is over the random phases $\bar{\theta}\left(\boldsymbol{{\rm k}}\right)$ which determine the stochastic function $\zeta\left({\bf r}\right)$.

The stochastic decompositions of the density matrix and of the Coulomb potential are then combined to give

[TABLE]

where the combined exchange-operator stochastic amplitude is simply $\xi({\bf r})=\eta({\bf r})\zeta({\bf r})$. The average is done now over $N_{\xi}$ random states; each sampling (labeled by $j$) of random $\xi({\bf r})$ amounts to a simultaneously choosing (independently) both the phases $\theta({\bf k})$ and the random vector $\bar{\eta}({\bf{\bf r})},$ so $\xi^{j}({\bf r})=\eta^{j}({\bf r})\zeta^{j}({\bf r})$.

Eq. (10) is formally exact if the number of states $N_{\xi}\to\infty$. For any finite number $N_{\xi}$ there is a statistical error proportional to $1/\sqrt{N_{\xi}}$, but since the long-range exchange vector is not numerically large this error is small even when $N_{\xi}$ is only a few hundreds. Further details, such as the supplementary use of a deterministic HOMO/LUMO when extracting $\gamma,$ are given in Ref. Neuhauser et al., 2015.

Note that in this present approach the only operator which is stochastically sampled is the long-range exchange. The density is still sampled deterministically from the eigenstates, $n({\bf r})=\sum_{i}f_{i}|\phi_{i}\left({\bf r}\right)|^{2}$, and the DFT cost is similar to that of traditional deterministic DFT for semilocal functionals. In practice, one can use any usual DFT algorithm to iteratively solves $H\phi_{i}=\varepsilon_{i}\phi_{i}$ for the occupied states, with $H$ constructed from $n({\bf r})$ and from $X_{{\rm nl}}^{\gamma}$ in Eq. (10).

The statistical errors in this mixed approach, where only the exchange is sampled stochastically, are much smaller than in our fully-stochastic DFT approachBaer, Neuhauser, and Rabani (2013); Neuhauser, Baer, and Rabani (2014) where the stochastic orbitals $\eta({\bf r})$ were also used to sample the local density (i.e., where we use $n({\bf r})=\left\{|\eta({\bf r})|^{2}\right\},$ or more generally $n({\bf r})=n_{0}({\bf r})+\left\{|\eta({\bf r})|^{2}\right\},$ where $n_{0}({\bf r})$ is a fragment density). In the fully-stochastic approach the eigenstates do not need to be determinedBaer and Neuhauser (2012); Baer, Neuhauser, and Rabani (2013); Neuhauser, Rabani, and Baer (2013) so it formally scales linearly with system size; here the scaling goal is more modest, just to reduce the cost to that of traditional DFT (and later TDDFT) with only local and semilocal potentials.

II.4 Stochastic TDDFT with non-local exchange

Following the DFT stage with stochastic exchange, we turn to the implementation of stochastic-exchange in TDDFT. There is no need to use the same exact methodology for the stochastic exchange in the TDDFT as in the DFT stage. Here, we follow Rabani, Baer, and Neuhauser (2015) and use a different sampling of the stochastic density matrix at each time step. Specifically, at each time-step we represent the density matrix as an average over stochastic vectors, where each one is constructed from the occupied eigenstates:

[TABLE]

where $\alpha_{i}(t)\in[0,2\pi]$ is a random phase. Thus, each $\beta$ is a stochastic vector created using a distinct set of random phases $\left\{\alpha_{i}\right\}$, and a different set of random phases is taken as each time step. Obviously $\left\{\beta({\bf r},t)\beta^{*}({\bf r^{\prime}},t)\right\}=\rho({\bf r},{\bf r^{\prime}},t).$

Note that the $\eta$ and $\beta$ vectors have a similar meaning; the former is used for the initial time-independent stage, the latter for TDDFT. We use a different symbol to emphasize that the number of such stochastic vectors is different in DFT and TDDFT. Specifically, since each time-step is small the effect of the stochastic exchange per time step is numerically small, so that it is sufficient to use only a small number ($N_{\beta}$) of vectors in each time step. For that reason, we have not done a stochastic resolution of the Coulomb kernel for the time-dependent exchange, which is formally calculated now as

[TABLE]

Thus the cost of evaluating Eq. 12 is only $N_{\beta}$-times more expensive than evaluating the Hartree term.

Since Eq. 12 is evaluated stochastically, $H\left(t\right)$ exhibits fluctuations even when there is no perturbation. For a linear response in $\Delta$, we therefore need to propagate two equations, with and without the perturbation:

[TABLE]

where $\phi_{i}^{\Delta}\left({\bf r},t=0\right)=\phi_{i}^{\Delta=0}\left({\bf r},t=0\right)=\phi_{i}\left({\bf r}\right)$. $H^{\Delta}\equiv H[n^{\Delta}(t),\rho^{\Delta}(t)]$ and $H^{\Delta=0}\equiv H^{\Delta}\equiv H[n^{\Delta=0}(t),\rho^{\Delta=0}(t)]$ have the same functional dependence on the density matrix, but since the time-dependent solutions of Eqs. 13 and 14 are different, we distinguish the Hamiltonians by superscripts. The time evolution of $\left\{\phi_{i}^{\Delta=0}\right\}$ stems purely from the stochastic fluctuations in $H^{\Delta=0}\left(t\right)$, as no external perturbing potential is applied. This fluctuation also induces time-dependence in the charge density $n^{\Delta=0}\left({\bf r},t\right)$ which needs to be subtracted when calculating the induced dipole

[TABLE]

from which the frequency-dependent absorption follows.

II.5 Numerical propagation of TDDFT with stochastic exchange

We use a split operator approach for the numerical propagation of the TDDFT equation with stochastic exchange. As usual, the perturbation is first applied at $t=0$ (and we again omit below the $\Delta$ superscript):

[TABLE]

and we then split the propagation of the non-local exchange and the remainder of the Hamiltonian,

[TABLE]

The short time kinetic+potential propagator (the non-$X_{{\it{\rm nl}}}^{\gamma}$ part in Eq. (17)) is itself calculated with a usual split operator evolution which will not be reviewed here, while $e^{-iX_{{\rm nl}}^{\gamma}\frac{dt}{2}}$ is evaluated extremely simply as

[TABLE]

where $N_{\phi}$ is a time-dependent normalization constant, i.e., $N_{\phi}^{-1}=||\left(1-iX_{{\rm nl}}^{\gamma}\frac{dt}{2}\right)|\phi\rangle||$. Since the normalization is dependent on the initial vector, Eq. (18) is slightly non-linear but this is of little practical consequence. The primitive approach of Eq. (18) is sufficient since the time-steps used are generally small, typically $dt=0.05$ a.u., i.e., around 1 as.

III Results

For all systems studied here, we first perform a ground-state DFT calculation and obtain the range-separation parameter $\gamma$ by enforcing the piecewise linearity condition for the total energy; this ensures that the HOMO is the same as the ionization energy.

III.1 Validation of the method using $\boldsymbol{{\rm PH_{3}}.}$

The smallest system studied is a PH3 molecule. Here, a deterministic DFT calculation was performed using experimental molecular structureHerzberg (1966) and the valence electronic states were computed with Troullier-Martins pseudopotentials.Troullier and Martins (1991) The total energy and the eigenvalues were converged to 5meV with a real space grid of $64\times 64\times 64$ points and a $0.4\,a_{0}$ grid spacing. Note that small molecular systems require, in general, a large range separation parameter and converge slower with the grid size and spacing compared to large systems. Through the tuning procedure,Stein et al. (2010); Kronik et al. (2012); Neuhauser et al. (2015) we found that $\gamma=0.37\,a_{0}^{-1}$; the resulting ionization potential (i.e., the negative of the HOMO energy) is $10.4$ eV, in excellent agreement with experiment ($10.6$ eV - Ref. Cowley et al., 1982).

The LUMO (obtained with the same range separation parameter) is barely bound, by slightly less than 0.1eV, but experimentally PH3 does not form a stable anion so the LUMO energy should be non-negative.

Using the optimally tuned BNL functional, the optical cross-section $\sigma\left(\omega\right)$ was obtained (see Eq. (A.6)) by deterministic and stochastic real time propagations, and the results are shown in Fig. 1. The TDDFT equations were propagated for a total time of 24 fs which provides a spectral resolution of $\sim$170 meV. The computed absorption cross-section has a first peak at $E_{1}=7.1$ eV, in excellent agreement with the first experimental peak at $7.0$ eV.

The exciton binding energy is defined as

[TABLE]

where $E_{g}$ is the fundamental band gap taken as the difference between HOMO and LUMO energies, i.e., $E_{g}=10.3$ eV. The predicted PH3 binding energy is thus $E_{b}=3.2$ eV, close to the experimental value of $3.4$ eV.111For PH3, we assume that the fundamental band gap coincides with the ionization potential as the molecule does not bind an extra electron

The overall absorption maximum is at $11.5$ eV, in good agreement with experiment ($12.0$ eV), though the latter exhibits large peak widths. At higher frequencies, the TDDFT spectrum has multiple local maxima (e.g., at $16.5$ and $20.4$ eV) that in the experiment only appear as shoulders. This is because experimental measurements cannot be resolved at energies higher than the ionization threshold ($10.4$ eV).

The stochastic decomposition of the time-dependent exchange (sTDDFT

• Eq. 12) reproduces the deterministic results already for $N_{\beta}=2$ (note that this is half the number used in deterministic exchange, which involves four valence states). We also checked that the spectrum does not change when the perturbation is varied in the range $\Delta=10^{-4}-10^{-3}$ a.u., as expected in linear response; this was also checked for the phosphorene sheets, discussed next.

III.2 Phosphorene Sheets: DFT

The major advantage of stochastic approaches is their applicability to large systems. We demonstrate this feature now on a set of 2D phosphorene sheets of increasing sizes, derived from a black phosphorus crystal structure.Cartz et al. (1979) The sheets were passivated with H atoms on the rims; each P atom that would have been bound to two H atoms was removed, resulting in a compact sheet geometry. Note that the two in-plane directions in the phosphorene sheet are traditionally labeled as armchair and zig-zag. With a kinetic energy cutoff of 26 $E_{h}$ and a real space grid with 208$\times$136$\times$40 points and a 0.6 $a_{0}$ spacing, the Kohn-Sham eigenvalues were converged to 10 meV.

The smallest sheet is $0.6\times 0.8$ nm and has 112 valence electrons. For its ground state DFT calculations we employed both the deterministic and stochastic formulation of the exchange operator $X_{{\rm nl}}^{\gamma}$. The optimally tuned range-separation parameter for this sheet is $\gamma=0.10\,a_{0}^{-1}$, and the stochastic eigenvalues converge slowly with the number of stochastic states, so $N_{\xi}\sim 1600$ is required to yield a statistical error of $<0.05$ eV.

In addition to the small sheet, we considered two larger sheets, $1.3\times 2.1$ nm and $3.1\times 4.3$ nm (labeled “medium” and “big”), with 478 and 1958 valence electrons respectively. For these larger sheets, the exchange operator was calculated purely stochastically as the deterministic calculation would have been very expensive. The range separation parameter gradually decreases with system size as in other 1D and 3D systems,Stein et al. (2010); Körzdörfer et al. (2011); Vlček et al. (2016) so $\gamma=0.09a_{0}^{-1}$ in the medium-size sheet and $0.05a_{0}^{-1}$ in the largest one. As the range separation parameter decreases with system size, the long-range exchange operator $X_{{\rm nl}}^{\gamma}$ is numerically smaller and its stochastic representation has therefore a small absolute statistical error. Hence, the largest system requires a smaller value of $N_{\xi}$ in the ground state calculations. Namely, the eigenvalues are converged to $<50$ meV with $N_{\xi}=1600$ stochastic states for the medium-size sheet and only $N_{\xi}=400$ for the largest one.

The stochastic-exchange DFT yields fundamental gaps $E_{g}$ that decrease with system size: $E_{g}=3.9$, $3.1$ and $1.7$ eV for the three sheets respectively. The large-sheet result is very similar to the HSE hybrid functional prediction of $1.5$ eV.Qiao et al. (2014) Comparing to with previous periodic $G_{0}W_{0}$ calculations (with a PBE starting point), we find that the largest sheet is in rough but not perfect agreement with the 2D periodic $G_{0}W_{0}$ fundamental gap of 2.08 eV.Li et al. (2017); Qiu, da Jornada, and Louie (2017) The difference from the $GW$ result could be due to the scalar (i.e., non-directional) nature of the range separated parameter, which ignores the difference between the effective interactions in the in- and out-of plane directions Vlček and Neuhauser , and perhaps also due to the approximate nature of the $G_{0}W_{0}$ itself.

III.3 Phosphorene sheets: TDDFT results

We next discuss the stability of TDDFT simulations for the sheets and the resulting optical spectra. Since the wavefunction is only incremented gradually, by $dt/2=0.025$ a.u, the statistical fluctuation introduced by each stochastic decomposition of the exchange operator is significantly smaller than for the ground state calculation. Hence, as mentioned, a small $N_{\beta}$ is sufficient for the time-dependent calculation, so the short-time results are fairly accurate already for $N_{\beta}=1-4$ . However, $N_{\beta}$ influences the total time of the simulations since due to statistical fluctuations the propagation eventually becomes unstable. We verified that this instability is not influenced by the time step, grid size and kinetic-energy cut-off.

We noticed the instability phenomena already in our original stochastic TDDFT approachGao et al. (2015) where, unlike here, we propagated only a few stochastic combination of eigenstates (i.e., several $\beta({\bf r},t))$, and constructed from them the density as $n({\bf r},t)\simeq N_{\beta}^{-1}\sum_{\beta}|\beta({\bf r},t)|^{2}$). That approach is extremely efficient for short time simulations (where the plasmon response of systems with thousands of electrons is accurately modeled by circa ten propagated states), but is limited to short times since the propagation eventually becomes unstable. Here, since all occupied eigenstates are propagated, and the density is constructed from all of them, the propagation is fairly stable for longer times.

Specifically, for the small and medium sheets, with 112 and 478 valence electrons, the propagation was carried up to $1000$ a.u ($\sim 24$ fs) without stability issues, using $N_{\beta}=$2 . This is in line with our previous simulations (Ref. Rabani, Baer, and Neuhauser, 2015) which used a damped exchange (reduced by 80%) in 3D and were stable with $N_{\beta}=1$. However, for the large sheet (with 1958 electrons) the time evolution became numerically unstable after $260$ a.u. ($\sim 7$ fs) even with $N_{\beta}$=6. The instability is for two reasons; first, the rapid oscillations of the density in the direction perpendicular to the phosphorene sheet. Unlike 3D systems, the response here is highly anisotropic and this appears to enhance the stochastic noise in the time propagation. Hence, for large 2D sheets the value of $N_{\beta}$ needs to be increased. Further, the form of Eq. (12) is oscillatory even for $\Delta=0$; in future publications we would use a less oscillatory form analogous to that in Ref. Rabani, Baer, and Neuhauser, 2015 and, in addition, would use a fully separable form of the TDDFT calculations, analogous to Eq. (10) for DFT with exchange.

Next, we turn to discuss the individual results. The time evolution of the induced dipole for the three phosphorene sheets is shown in Fig. 2. The absorption cross sections per-electron are shown in Fig. 3 for the zig-zag and armchair directions. The spectra are strongly anisotropic, in agreement with experimental data and $GW$/BSE calculations on 2D periodic sheets.Liu et al. (2014); Xia, Wang, and Jia (2014); Wang et al. (2015); Li et al. (2017) With rising system size there are diminished quantum confinement effects so the fundamental band gap decreases and the absorption maximum therefore gradually shifts to lower energies. As the number of valence electrons increases the spectra also smoothens.

On the right panel of Fig. 3 we zoom on the absorption spectrum below the ionization potential (which only slowly decreases with system size, from 5.9 to 5.4 eV). The absorption cross section decreases rapidly at lower frequencies, but several local maxima are evident. Some of these local maxima are below the band gap energy ($E_{g}$) and therefore correspond to bound electron-hole pairs. These excitonic peaks appear only for the armchair direction (due to the strong anisotropy of the optical response); this feature was seen in previous calculations for periodic phosphorene and was also seen experimentally.Liu et al. (2014); Tran et al. (2014); Li et al. (2017)

Excitonic peaks are usually sharp and have a high intensity, indicating long-lived quasiparticle states. The maxima in Fig. 3 are however broadened due to the finite simulation time (24 fs for the two small systems and 7 fs for the largest one). For the small sheet, the excitonic peaks are well-separated but have relatively low intensity. The position of the first absorption peak maximum ($E_{1}$) changes with increasing system size from $1.8$ to $1.6$ eV. The latter is in good agreement with the experimental value and $GW$/BSE estimates, $1.7$ and $1.6$ eV, respectively Li et al. (2017); Qiu, da Jornada, and Louie (2017) for a bulk 2D system.

The exciton binding energy (Eq. 19) decreases rapidly with system size from $2.1$ to $0.1$ eV. The strongest excitonic response (the largest amplitude of the $E_{1}$ peak) is found in the medium sized system, which also has a high exciton binding energy $E_{b}=1.65$ eV. As mentioned, however, the fundamental gap $E_{g}$ in the stochastic-exchange DFT is underestimated relative to $G_{0}W_{0}$ calculations. Therefore, the exciton binding energy for the largest sheet (0.1 eV) is much lower than predicted by $GW$/BSE calculations which give $E_{b}=0.48$ eV.Qiu, da Jornada, and Louie (2017)

Interestingly, when the phosphorene is encapsulated in dielectric media, the $GW$/BSE binding energy becomes small, $0.14$ eV Qiu, da Jornada, and Louie (2017), comparable to our TDDFT estimates of pristine (non-encapsulated) phosphorenes. The encapsulation causes strong screening above and below the 2D system (i.e., in the out-of-plane direction) leading to a big change in the $G_{0}W_{0}$ gap $E_{g}$ (from $2.08$ to $1.62$ eV, similar to our pristine DFT gap), while the position of the first excitonic peak, $E_{1}$ remains practically unaffected.Qiu, da Jornada, and Louie (2017) The difference between our results and experiment and $GW$/BSE points to a problem in describing 2D materials with range-separated potentials, Vlček and Neuhauser which should in principle account for the anisotropy of the electron-electron interaction between the in- and out-of-plane directions. Until such a non-isotropic interaction is implemented in DFT and TDDFT, stochastic TDDFT with exchange can only be trusted as far as the exciton frequency, but calculations of the exciton binding would require a more accurate method than GKS-DFT for the quasiparticle gap.

IV Conclusions

In summary, we developed an efficient real-time TDDFT approach with stochastic long-range non-local exchange. The stochastic treatment decomposes the density matrix in TDDFT to an average over a product of random vectors $\beta$ in the space spanned by the occupied orbitals. It significantly reduces the computational cost as only a few stochastic states are needed at each time step. Further, the number of stochastic states varies only a little with the system size. Calculations for very large systems thus become feasible.

The resulting TDDFT with long-range non-local exchange includes the attractive electron-hole interaction that gives rise to exciton formation. Indeed, our TDDFT yields optical spectra that are in excellent agreement with experiment. For small systems, where deterministic calculations are affordable, the stochastic and deterministic results agree.

We demonstrated that our method is applicable for extremely big systems using a set of phosphorene sheets containing up to $\sim$2000 valence electrons. The largest system was compared to experiments and previous calculations on infinite phosphorene sheets; the analysis confirms that the range-separated hybrid functional successfully predicts optical spectra even with strong excitonic signatures.

V Acknowledgments

D.N. acknowledges support by NSF grant CHE-1763176. RB acknoledges the US-Israel Binational Fund grant no. BSF2015687. The authors would like to acknowledge helpful discussions with Eran Rabani. The calculations were performed as part of the XSEDETowns et al. (2014) computational project TG-CHE180051.

Appendix I: Photoabsorption cross section

Here we overview for completeness the well-known expression of the photoabsorption cross-section as a Fourier transform of a real-time dipole correlation function.

The absorption cross section, $\sigma(\omega)$, is given in linear response as Onida, Reining, and Rubio (2002):

[TABLE]

where tilde is used occasionally to denote quantities in frequency domain, $\delta\tilde{v}\left({\bf r},\omega\right)$ is dynamical external potential and $\chi$ is the electronic reducible polarizability, which is given in the time domain as:

[TABLE]

where $\delta n\left({\bf r},t\right)$ is the induced charge density at a point ${\bf r}$ and time $t$. The response function is causal, i.e., $t>t^{\prime}$. For absorption of polarized light along a unit vector $\hat{{\rm\boldsymbol{e}}}$ we apply $\delta v$ as a dipole potential. The cross section is then Onida, Reining, and Rubio (2002)

[TABLE]

Here, $\sigma$ is calculated from real-time linear-response. Specifically, the first stage is to apply a dipole external potential perturbation

[TABLE]

where $\Delta$ is the perturbation strength, and an instantaneous perturbation is applied at $t=0$ allowing to probe the response at all frequencies. The perturbation potential is applied to all occupied eigenstates, which are then propagated in time. The resulting oscillations of the induced charge density (Eq. A.2) are then used to find the dipole auto-correlation,

[TABLE]

The absorption cross section is finally a Fourier transform of the dipole auto correlation:

[TABLE]