Basis-Independent Spectral Methods for Non-linear Optical Response in Arbitrary Tight-binding Models
NNon-linear optical response of non-interacting systems with spectral methods
S. M. João
Centro de Física das Universidades do Minho e Porto andDepartamento de Física e Astronomia, Faculdade de Ciências,Universidade do Porto, 4169-007 Porto, Portugal
J. M. V. P. Lopes
Centro de Física das Universidades do Minho e PortoDepartamento de Engenharia Física, Faculdade de Engenharia andDepartamento de Física e Astronomia, Faculdade de Ciências,Universidade do Porto, 4169-007 Porto, Portugal
A perturbation expansion is developed using the Keldysh formalism to obtain the nonlinear opticalconductivity of non-interacting systems with special focus on tight-binding models. This basis-independent approach provides the objects required by numerical spectral methods such as theKernel Polynomial Method (KPM), allowing for a straightforward numerical implementation of thelinear and nonlinear optical conductivity. Using the real-space basis, we present results for thesecond-order conductivity in hexagonal Boron Nitride with vacancies and with Anderson disorder.
I. INTRODUCTION
Since the advent of the laser in 1960, the field of non-linear optics has received considerable interest. The fol-lowing year marked the beginning of a systematic studyof this field, as in 1961 P. Franken was able to demon-strate second harmonic generation (SHG) [7] experimen-tally. This opened the gateway to a whole new plethoraof phenomena, such as optical rectification (1962) [2] andhigher harmonic generation (1967) [15]. The laser wasable to provide the powerful electric fields needed to seethese nonlinear phenomena. The optical properties ofcrystals have been studied extensively throughout thelast 30 years [21] and have recently received consider-able interest due to the strong non-linear properties oflayered materials such as Graphene [3, 9, 14].Several approaches have been developed to obtain thenon-linear response up to arbitrary order of a crystallinesystem subject to an external field. Among those, wemay find generalizations of Kubo’s formula for higher or-ders and perturbation expansions for the density matrixin both the length gauge and the velocity gauge [21, 22].Due to the complexity of the problem, the bulk of thatwork was done using translation-invariant systems. De-spite covering a broad range of systems, this approachis limited in its scope, as it cannot deal with impurities,border effects, or even magnetic fields within the usualPeierls framework [4]. In order to simulate more realis-tic systems, this restriction has to be lifted. The goal isthen to re-express the same formulas for the linear andnonlinear optical conductivities in a framework that al-lows real-space calculations. In this work, we extend thework of Passos et al [17] implementing the velocity gaugemethodology within the Keldysh formalism [11]. Thisway, we develop a general perturbation procedure to dealwith non-interacting fermion systems at finite tempera-ture coupled to a time-dependent external field. Through careful categorization of all these contributions, we pro-vide a systematic procedure to find the objects neededto calculate the conductivity at any order. These ob-jects are expressed with no reference to a specific basis.The critical point here is that the mathematical objectsprovided by our perturbation expansion are precisely theones required by the numerical spectral methods we use.This fact, combined with our diagrammatic approach,provides a straightforward way to implement the numer-ical calculation of the nonlinear optical conductivity fora wide range of materials. This methodology is imple-mented in KITE, an open-source software developed byourselves [13] to calculate the nonlinear optical propertiesof a very broad range of tight-binding models.This paper is structured as follows. In Section II, weintroduce the Keldysh formalism to show how the currentwill be calculated and define the fundamental mathemat-ical objects of our calculations. Section III applies theperturbation procedure to a tight-binding model. Theelectromagnetic field is added through the Peierls Sub-stitution. Then, a diagrammatic procedure is developedin order to deal with the large number of non-trivial con-tributions to the conductivity in each order. Section IVexplains how to use a Chebyshev expansion of the op-erators in the previous expression in order to be usedwith spectral methods. This provides the full formulathat may be directly implemented numerically. Finally,in Section V we apply the formalism developed in theprevious sections to several different systems. We startby comparing with known results obtained by k -space in-tegration [17, 22] of the nonlinear optical conductivity ofhexagonal Boron Nitride. Then, we add some disorder tothe system and show that our method gives the expectedresult. Finally, the numerical efficiency and convergenceof our method is assessed. a r X i v : . [ c ond - m a t . o t h e r] M a r II. KELDYSH FORMALISM
The Keldysh formalism [6, 12, 20] is a general pertur-bation scheme describing the quantum mechanical timeevolution of non-equilibrium interacting systems at finitetemperature. It provides a concise diagrammatic repre-sentation of the average values of quantum operators.This formalism does not rely on any particular basis,which is a critical feature for this paper. The math-ematical objects obtained in our expansion are in theform required by numerical spectral methods, providinga straightforward formula for the numerical calculationof nonlinear optical response functions. In this section wewill introduce the definitions of the objects used through-out the paper and show how to expand the Green’s func-tions for fermions [11] with this formalism.
A. Definitions
1. Green’s functions
To use the Keldysh formalism for fermions, we needthe definitions of the time-ordered, lesser, greater andanti-time-ordered Green’s functions. Respectively, iG Tab ( t, t (cid:48) ) = (cid:68) T (cid:104) c a ( t ) c † b ( t (cid:48) ) (cid:105)(cid:69) (1) iG
2. Expected value of an operator
The expected value of the current J ( t ) (or any one-particle operator) may be evaluated with resort to these Figure 1: Diagrammatic representation of the expected valueof the current operator in Fourier space. The horizontalstraight line ending in a circle is the lesser Green’s functionand the wavy line beginning in a circle represents the currentoperator.
Green’s functions by tracing over its product with theperturbed lesser Green’s function: J ( t ) = (cid:68) ˆ J ( t ) (cid:69) = − Tr (cid:104) ˆ J ( t ) iG < ( t, t ) (cid:105) . (7)The Fourier transform [24] of J ( t ) is shown diagram-matically in Figure 1. The circles stand for the full, per-turbed operators in the presence of an external field.
3. Conductivity
We use the same definition for the nonlinear opticalconductivity as in [17, 22]: J α ( ω ) = σ αβ ( ω ) E β ( ω ) + ˆ d ω π ˆ d ω π × (8) σ αβγ ( ω , ω ) E β ( ω ) E γ ( ω ) 2 πδ ( ω + ω − ω ) + · · · where E α is the component of the electric field along the α direction and the repeated indices are assumed to besummed over. The coefficients of this expansion are theconductivities at each order in the expansion. The nextsection is devoted to finding the perturbation expansionof G < . In this paper we are dealing with tight-bindingmodels, in which case the current operator will itself bea power series of the external field. B. Non-interacting electronic systems
Our system is described by the many-particle time-dependent Hamiltonian H ( t ) = H + H ext ( t ) . (9)where H is an Hamiltonian that we can solve exactlyand H ext ( t ) is the time-dependent external perturbation.Here we restrict ourselves to non-interacting Hamiltoni-ans since we’re dealing with non-interacting electrons.These operators are expressed in terms of their single-particle counterparts as H ext ( t ) = (cid:88) ab [ H ext ( t )] ab c † a ( t ) c b ( t ) (10) H = (cid:88) ab [ H ] ab c † a ( t ) c b ( t ) . (11)The expansion of the perturbed lesser Green’s function G < will be expressed in terms of the unperturbed Green’sfunctions g > , g < , g R and g A in Fourier space: i ˜ g < ( ω ) = − πf ( (cid:126) ω ) δ ( ω − H / (cid:126) ) (12) i ˜ g > ( ω ) = 2 π [1 − f ( (cid:126) ω )] δ ( ω − H / (cid:126) ) (13) i ˜ g R ( ω ) = iω − H / (cid:126) + i + (14) i ˜ g A ( ω ) = iω − H / (cid:126) − i + , (15)where f ( (cid:15) ) = (cid:0) e β ( (cid:15) − µ ) (cid:1) − is the Fermi-Dirac distri-bution, β is the inverse temperature and µ is the chemicalpotential. The Keldysh formalism and Langreth’s rulesprovide the perturbation expansion of G < [11]. Defin-ing V ( t ) = ( i (cid:126) ) − H ext ( t ) , the zeroth-order term in theexpansion is i ˜ G < (0) ( ω ) = ˆ d ω i ˜ g < ( ω ) δ ( ω ) (16)and the first-order one is i ˜ G < (1) ( ω ) = ˆ d ω (2 π ) (2 π ) δ ( ω − ω − ω ) (17) × δ ( ω + ω − ω ) (cid:104) i ˜ g R ( ω ) ˜ V ( ω ) i ˜ g < ( ω )++ i ˜ g < ( ω ) ˜ V ( ω ) i ˜ g A ( ω ) (cid:105) . ´ d n ω ··· n is a shorthand for ´ · · · ´ d ω · · · d ω n . Thesecond-order term is i ˜ G < (2) ( ω ) = ˆ d ω ··· (2 π ) (2 π ) δ ( ω + ω − ω ) × δ ( ω − ω − ω ) δ ( ω − ω − ω ) × (cid:104) i ˜ g R ( ω ) ˜ V ( ω ) i ˜ g R ( ω ) ˜ V ( ω ) i ˜ g < ( ω )+ i ˜ g R ( ω ) ˜ V ( ω ) i ˜ g < ( ω ) ˜ V ( ω ) i ˜ g A ( ω ) (18) + i ˜ g < ( ω ) ˜ V ( ω ) i ˜ g A ( ω ) ˜ V ( ω ) i ˜ g A ( ω ) (cid:105) . Diagrammatically, the expansion of iG < ( ω ) is repre-sented by Figure 2. Each wavy line ending in a circle rep-resents an external perturbation ˜ V . There are three dif-ferent types of Green’s functions that may appear in theseexpansions, with a certain regularity: a lesser Green’sfunction ˜ g < , which is always present, retarded Green’s = +++ + + Figure 2: Diagrammatic representation of the lesser Green’sfunction. functions ˜ g R and advanced Green’s functions ˜ g A . Dia-grammatically, ˜ g < is represented by a dashed line whilethe solid lines represent retarded or advanced Green’sfunctions. To identify whether a line represents a re-tarded or advanced Green’s function, one needs to readthe diagram and identify the position of the lesser Green’sfunction and the outgoing line. Reading clockwise (anti-clockwise) until finding the outgoing line, there can onlybe advanced (retarded) Green’s functions. In each in-tersection, the corresponding external perturbation ˜ V isinserted. An exception is made for the intersection withthe line representing ω , as it still needs to be contracted.If the external perturbation were a simple external field E ( t ) , then the coupling would be H ext ( t ) = e E ( t ) · r andthe previous expressions coupled with eq. (7) would suf-fice. Now we will turn to tight-binding Hamiltonians, forwhich the external coupling is actually an infinite seriesof operators due to the way the electromagnetic field isintroduced. This affects not only the V operators butalso the expression for the current operator. III. TIGHT-BINDING HAMILTONIAN WITHEXTERNAL ELECTRIC FIELD
Tight-binding models provide a simple framework withwhich to calculate transport quantities. This frameworkcan be used to express structural disorder in the system,whilst Peierls’ substitution [18] adds an electromagneticfield as an external perturbation. Despite the simplicityof this procedure, the addition of an electromagnetic fieldthrough a phase factor yields an infinite series of H ext .In this section, we obtain the expression for H ext andshow how the expansions of the previous sections may beused to obtain the nonlinear optical conductivity. This isentirely analogous to the way the external perturbation isintroduced with the velocity gauge in the work of Passoset al [17]. A. Series expansion
Let’s consider the following tight-binding Hamiltonian: H = (cid:88) R i , R j (cid:88) σ ,σ t σ σ ( R i , R j ) c † σ ( R i ) c σ ( R j ) . (19)The R i represent the lattice sites and the σ i the otherdegrees of freedom unrelated to the position, such as theorbitals and spin. The electromagnetic field is introducedthrough Peierls’ substitution: t σ σ ( R i , R j ) → e − ie (cid:126) ´ R i R j A ( r (cid:48) ,t ) · d r (cid:48) t σ σ ( R i , R j ) . (20)To introduce both a static magnetic field and a uniformelectric field, we use the following vector potential: A ( r , t ) = A ( r ) + A ( t ) . (21)The electric and magnetic fields are obtained from E ( t ) = − ∂ t A ( t ) and B ( r ) = ∇ × A ( r ) . The introduc-tion of the magnetic field only changes the t σ σ ( R i , R j ) without introducing a time dependency. Therefore, wemay assume that a magnetic field is always present with-out any loss of generality for the following discussionwhile keeping in mind that its introduction broke trans-lation invariance. Since the magnetic field only affectsthe hopping parameters, from now on, the term in thevector potential that provides the electric field will bedenoted by A ( t ) . The external perturbation is obtainedby expanding the exponential in eq. 20 and identifyingthe original Hamiltonian.
1. Expansion of the external perturbation
Expanding the exponential in eq. 20 yields an infiniteseries of operators for the full Hamiltonian: H A ( t ) = H + H ext ( t ) (22)from which we identify, after a Fourier transform, ˜ V ( ω ) = ei (cid:126) h α ˜ A α ( ω ) + e i (cid:126) h αβ ˆ d ω (cid:48) π ˆ d ω (cid:48)(cid:48) π × (23) ˜ A α ( ω (cid:48) ) ˜ A β ( ω (cid:48)(cid:48) ) 2 πδ ( ω (cid:48) + ω (cid:48)(cid:48) − ω ) + · · · . Repeated spatial indices are understood to be summedover. We have defined ˆ h α ··· α n = 1( i (cid:126) ) n [ˆ r α , [ · · · [ˆ r α n , H ]]] (24) = + + + ... Figure 3: Diagrammatic representation of the external per-turbation. = + + + ...
Figure 4: Diagrammatic representation of the current oper-ator. The single small circle is to be understood as a Diracdelta. where ˆ r is the position operator. In first order, ˆ h α isjust the single-particle velocity operator. In Figure 3, wesee how the diagrammatic representation of the externalperturbation unfolds into an infinite series of externalfields. The wavy line represents ˜ A and the number ofexternal fields connected to the same point is the numberof commutators in eq. 24.
2. Expansion of the current
The current operator is calculated directly from theHamiltonian, using ˆ J α = − Ω − ∂H/∂A α ( Ω is the volumeof the sample), which also follows a series expansion dueto the presence of an infinite number of A ( t ) in H ext : ˆ J α ( t ) = − e Ω (cid:16) ˆ h α + e ˆ h αβ A β ( t ) ++ e
2! ˆ h αβγ A β ( t ) A γ ( t ) + · · · (cid:19) . (25)Figure 4 depicts the diagrammatic representation ofthis operator in Fourier space.The complexity of this expansion becomes clear. Ineq. 7, both the current operator and the Green’s func-tions follow a perturbation expansion. Furthermore, eachinteraction operator in every one of the terms in theGreen’s function expansion also follows a similar expan-sion. We now have all the objects needed for the pertur-bative expansion of the conductivity. B. Perturbative expansion of the conductivity
In the previous sections we laid out the expressionsfor each individual operator in our expansion and repre-sented their corresponding diagrammatic depictions. Inthis subsection, we put together all the elements of theprevious sections to provide the full diagrammatic rep-resentation of the first and second-order conductivities. + ++++ + ++ + a)b)
Figure 5: Expansion of the expected value of the conductivityin (a) first and (b) second order.
This expansion closely resembles that of [16] but has sev-eral differences due to the usage of these specific Green’sfunctions. The only thing left to do is to replace theperturbed objects in the diagrammatic representation ofthe expected value of the current operator by their ex-pansions. It is straightforward to see how the diagramsfit together in Figure 5, which shows all the contributing diagrams up to second order.Obtaining the conductivity from the current is a mat-ter of expressing the frequencies ω (cid:48) , ω (cid:48)(cid:48) and ω (cid:48)(cid:48)(cid:48) in termsof ω , ω and ω and using E ( ω ) = iω A ( ω ) . The Diracdelta in eq. 8 simply means that ω is to be replacedby the sum of external frequencies entering the diagram.Thus, the n-th order conductivity may be found usingthe following rules:1. Draw all the diagrams with n wavy lines comingin the diagram, one going out and one dashed in-terconnecting line. Integrate over the internal fre-quencies and ignore the conservation of momentumin the vertex containing ω , as that is already takeninto account by the Dirac delta in the definition ofthe conductivity.2. Reading clockwise starting from the vertex contain-ing ω , insert, by order, a generalized velocity opera-tor h α ··· α k at each vertex and a Green’s function ateach edge. Each α i is the label of a frequency lineconnecting to the vertex. If the edge is a dashedline, the Green’s function is ig < . All the edges be-fore that correspond to ig R and the ones after it to ig A . Trace over the resulting operator.3. Multiply by Ω − e n +1 (cid:81) nk =1 ( iω k ) − ( i (cid:126) ) − N , where n is the number of dashed lines and N is the numberof interconnecting lines. For each vertex, divide bythe factorial of the number of outgoing lines.Following these rules and replacing ig < by eq. 12, thefirst-order conductivity is found: σ αβ ( ω ) = ie Ω ω ˆ ∞−∞ d(cid:15)f ( (cid:15) ) Tr (cid:20) ˆ h αβ δ ( (cid:15) − H ) + 1 (cid:126) ˆ h α g R ( (cid:15)/ (cid:126) + ω ) ˆ h β δ ( (cid:15) − H ) + 1 (cid:126) ˆ h α δ ( (cid:15) − H ) ˆ h β g A ( (cid:15)/ (cid:126) − ω ) (cid:21) . (26)Similarly, for the second-order conductivity: σ αβγ ( ω , ω ) = 1Ω e ω ω ˆ ∞−∞ d(cid:15)f ( (cid:15) ) Tr (cid:20)
12 ˆ h αβγ δ ( (cid:15) − H ) + 1 (cid:126) ˆ h αβ g R ( (cid:15)/ (cid:126) + ω ) ˆ h γ δ ( (cid:15) − H ) (27) + 1 (cid:126) ˆ h αβ δ ( (cid:15) − H )ˆ h γ g A ( (cid:15)/ (cid:126) − ω ) + 12 (cid:126) ˆ h α g R ( (cid:15)/ (cid:126) + ω + ω ) ˆ h βγ δ ( (cid:15) − H ) + 12 (cid:126) ˆ h α δ ( (cid:15) − H )ˆ h βγ g A ( (cid:15)/ (cid:126) − ω − ω )+ 1 (cid:126) ˆ h α g R ( (cid:15)/ (cid:126) + ω + ω ) ˆ h β g R ( (cid:15)/ (cid:126) + ω ) ˆ h γ δ ( (cid:15) − H ) + 1 (cid:126) ˆ h α g R ( (cid:15)/ (cid:126) + ω ) ˆ h β δ ( (cid:15) − H )ˆ h γ g A ( (cid:15)/ (cid:126) − ω )+ 1 (cid:126) ˆ h α δ ( (cid:15) − H )ˆ h β g A ( (cid:15)/ (cid:126) − ω ) ˆ h γ g A ( (cid:15)/ (cid:126) − ω − ω ) (cid:21) . The procedure is exactly the same for the n -th order conductivity, which will have n − ( n + 2) diagrams. Thehigher-order expansions will not be obtained because arealistic computation of physical quantities with thoseformulas using spectral methods would require tremen-dous computational power. This point will be furtherexplained in the next section. Despite these limitations,we have provided the required framework for obtainingthose formulas, which might be useful in the future. IV. SPECTRAL METHODS
From the previous section, it becomes clear that theonly objects needed to calculate the conductivity up toany order are the retarded and advanced Green’s func-tions, Dirac deltas and the generalized velocity operators.The meaning of these functions is obvious in the energyeigenbasis of H , but not in the position basis that weultimately want to use. Therefore, these objects are ex-panded in a truncated series of Chebyshev polynomials.This choice of polynomials allows for a very efficient andnumerically stable method of computing these functions[23]. The fact that the Dirac deltas and Green’s functionshave singularities means that their expansions will beplagued by Gibbs oscillations. Some methods like KPMadd a weight function [8, 19] to each term in the ex-pansion, effectively damping the oscillations. Let’s takeas an example the Green’s function. Using a kernel, asmore polynomials are added to the expansion, the expan-sion becomes closer to the exact Green’s function, and somore singular. Although the behavior approaches thatof a Green’s function as more polynomials are added, wecannot speak of convergence in the usual sense. In orderto assess the convergence properties of our method, wewill not use a kernel, but instead use a finite imaginarybroadening parameter inside the Green’s function, thatis i + → iλ . For a finite λ , the function is no longersingular and so we can expect the expansion to convergewithin a given accuracy after enough polynomials havebeen added. In this paper, we will use an exact decompo-sition of the Green’s function [5] in terms of Chebyshevpolynomials in order to be able to evaluate the conver-gence of our method. The term (cid:126) /λ may also be in-terpreted as a phenomenological relaxation time due toinelastic scattering processes and therefore may be ad-justed to reflect this fact.In this section, we introduce the Chebyshev polyno-mials and use them to expand the Green’s functions andDirac deltas of the previous formulas. This will cast thoseformulas into a more useful form. A. Expansion in Chebyshev polynomials
The Chebyshev polynomials of the first kind are a setof orthogonal polynomials defined in the range [ − , by T n ( x ) = cos ( n arccos ( x )) . (28)They satisfy a recursion relation T ( x ) = 1 (29) T ( x ) = x (30) T n +1 ( x ) = 2 xT n ( x ) − T n − ( x ) (31)and the following orthogonality relation ˆ − T n ( x ) T m ( x ) d x √ − x = δ nm δ n . (32)These polynomials may be used to expand functions ofthe Hamiltonian provided its spectrum has been scaledby a factor ∆ E to fit in the range [ − , . The onlyfunctions of the Hamiltonian appearing in the expansionare Dirac deltas and Green’s functions, which have thefollowing expansions in terms of Chebyshev polynomials: δ ( (cid:15) − H ) = ∞ (cid:88) n =0 ∆ n ( (cid:15) ) T n ( H )1 + δ n (33) g σ,λ ( (cid:15), H ) = (cid:126) (cid:15) − H + iσλ = (cid:126) ∞ (cid:88) n =0 g σ,λn ( (cid:15) ) T n ( H )1 + δ n (34)where ∆ n ( (cid:15) ) = 2 T n ( (cid:15) ) π √ − (cid:15) (35)and g σ,λn ( (cid:15) ) = − σi e − niσ arccos( (cid:15) + iσλ ) (cid:113) − ( (cid:15) + iσλ ) . (36)The function g σ,λ encompasses both retarded and ad-vanced Green’s functions in the limit λ → + : g + , + isthe retarded Green’s function and g − , + the advancedone. The operator part has been completely separatedfrom its other arguments. All the Dirac deltas andGreen’s functions may therefore be separated into a poly-nomial of H and a coefficient which encapsulates thefrequency and energy parameters. The trace in the con-ductivity now becomes a trace over a product of polyno-mials and ˆ h operators, which can be encapsulated in anew object, the Γ matrix: Γ α , ··· , α m n ··· n m = Tr N (cid:20) ˜ h α T n ( H )1 + δ n · · · ˜ h α m T n m ( H )1 + δ n m (cid:21) . (37)The upper indices in bold stand for any number of in-dices: α = α α · · · α N . Here we have used ˜ h α =( i (cid:126) ) N ˆ h α rather than ˆ h to avoid using complex num-bers when the Hamiltonian matrix is purely real in ournumerical simulations. It’s very important to keep inmind that these new operators are no longer hermitian.The commas in Γ separate the various ˜ h operators. N isthe number of unit cells in the sample being studied andensures that Γ is an intensive quantity. Some examples: Γ α,βγnm = Tr N (cid:20) ˜ h α T n ( H )1 + δ n ˜ h βγ T m ( H )1 + δ m (cid:21) (38) Γ αβn = Tr N (cid:20) ˜ h αβ T n ( H )1 + δ n (cid:21) (39) Γ α,β,γnmp = Tr N (cid:20) ˜ h α T n ( H )1 + δ n ˜ h β T m ( H )1 + δ m ˜ h γ T p ( H )1 + δ p (cid:21) . (40)The Γ matrix only depends on the physical system it-self as it is merely a function of the Hamiltonian and the ˜ h operators. The coefficients of the Chebyshev expan-sion may similarly be aggregated into a matrix, whichwe denote by Λ . Some examples: Λ n = ˆ ∞−∞ d(cid:15)f ( (cid:15) )∆ n ( (cid:15) ) (41) Λ nm ( ω ) = (cid:126) ˆ ∞−∞ d(cid:15)f ( (cid:15) ) (cid:2) g Rn ( (cid:15)/ (cid:126) + ω ) ∆ m ( (cid:15) )+∆ n ( (cid:15) ) g Am ( (cid:15)/ (cid:126) − ω ) (cid:3) (42)and Λ nmp ( ω , ω ) = (43) (cid:126) ˆ ∞−∞ d(cid:15)f ( (cid:15) ) (cid:2) g Rn ( (cid:15)/ (cid:126) + ω + ω ) g Rm ( (cid:15)/ (cid:126) + ω ) ∆ p ( (cid:15) )+ g Rn ( (cid:15)/ (cid:126) + ω ) ∆ m ( (cid:15) ) g Ap ( (cid:15)/ (cid:126) − ω )+∆ n ( (cid:15) ) g Am ( (cid:15)/ (cid:126) − ω ) g Ap ( (cid:15)/ (cid:126) − ω − ω ) (cid:3) In terms of these new objects, the conductivities be-come σ αβ ( ω ) = − ie Ω c (cid:126) ω (cid:34)(cid:88) n Γ αβn Λ n + (cid:88) nm Λ nm ( ω ) Γ α,βnm (cid:35) (44)in first order and σ αβγ ( ω , ω ) = ie Ω c ω ω (cid:126) (cid:34) (cid:88) n Λ n Γ αβγn + (cid:88) nm Λ nm ( ω ) Γ αβ,γnm + 12 (cid:88) nm Λ nm ( ω + ω ) Γ α,βγnm + (cid:88) nmp Λ nmp ( ω , ω ) Γ α,β,γnmp (cid:35) (45) in second order. Ω c is the volume of the unit cell. B. Considerations on the numerical storage of Γ Naturally, one cannot expect to sum the entire Cheby-shev series, so it has to be truncated at a certain numberof polynomials N max . Each of the entries in a Γ matrixrepresents a complex number. Numerically, this is rep-resented as two double-precision floating-point numbers,each taking up bytes of storage. The amount of stor-age needed to store a Γ matrix of dimension n is N n max .The number of Chebyshev polynomials needed to obtaina decent resolution depends heavily on the problem athand, but a typical number may be N max = 1024 . A one-dimensional Γ matrix would take up Kb of storage, atwo-dimensional matrix Mb and a three-dimensionalmatrix Gb. Three-dimensional matrices appear in thesecond-order conductivity. The third-order conductivitywould require a four-dimensional matrix and as such, Tb of storage. Numbers like these make it unrealistic togo beyond second order conductivity.
V. NUMERICAL RESULTS
In this section we showcase several examples, of in-creasing complexity, to compare our formalism with theliterature. Starting with Graphene, we compute the lin-ear optical conductivity and verify that it agrees per-fectly with the k -space formalism. Breaking the sublat-tice symmetry with hexagonal Boron Nitride (h-BN), weare able to obtain the photoconductivity and check thatit too agrees perfectly. This proves that our method isable to accurately reproduce the existing results. Then,we show two examples that cannot be reproduced easilywith the k -space formalism: second harmonic generation(SHG) in h-BN with Anderson disorder and vacanciesof varying concentration. Finally, the convergence prop-erties are evaluated and the efficiency of the method isdiscussed. A. Linear optical response in Graphene
Let a be the distance between consecutive atoms in thehoneycomb lattice. Then, the primitive vectors betweenunit cells are (see Figure 6) a = a (cid:16) √ , (cid:17) a = a (cid:32) √ , (cid:33) and the distance vectors between nearest neighbors are a a 𝛿 𝛿 𝛿 B A xy Figure 6: Honeycomb lattice and choice of primitive vectors. δ = a (cid:16) √ , − (cid:17) δ = a (0 , δ = a (cid:16) −√ , − (cid:17) . The area of the unit cell is Ω c = √ a . Start-ing from eq. 19, the Graphene Hamiltonian is ob-tained by invoking translational invariance of the unitcell t µν ( R m , R n ) = t µν ( R m − R n ) and t AB ( δ ) = t AB ( δ ) = t AB ( δ ) = − t. The remaining non-zero hopping integrals are foundby using t AB = t BA . The on-site energies t AA ( ) and t BB ( ) are taken to be zero. A factor of two is includeddue to spin degeneracy.These parameters were used to obtain the first-orderoptical conductivity for Graphene, as seen in Figure 7.We used a lattice with unit cells in each directionand Chebyshev moments in the expansion. The re-sulting plot is compared to the results obtained in [22]through k -space integration of a translation-invariantsystem. The curves are indistinguishable. B. Hexagonal Boron Nitride
The only difference relative to Graphene is found inthe on-site energies. Let t AA ( ) = ∆ / and t BB ( ) = − ∆ / . In this model, the one- and three-dimensionaland Γ matrices are identically zero[25], so the second-order conductivity may be calculated resorting only totwo-dimensional Γ matrices. The calculation is thus sim-plified tremendously because the conductivity reduces to σ αβγ ( ω , ω ) = ie Ω c ω ω (cid:126) × (cid:34)(cid:88) nm Λ nm ( ω ) Γ αβ,γnm + 12 (cid:88) nm Λ nm ( ω + ω ) Γ α,βγnm (cid:35) . Figure 7: First-order longitudinal yy conductivity forGraphene in units of σ = e / (cid:126) . Hopping parameter: t = 2 . eV, temperature: T = 2 . mK, chemical potential: µ = 0 . eV, broadening parameter: λ = 38 . meV, num-ber of Chebyshev moments used: M = 2048 , lattice size: L = 4096 × . The solid curves represent the optical con-ductivity obtained by KITE (real part in green, imaginary inblue). The superimposed dashed lines are obtained in [22]. The indices n, m are understood to be summed over.The photoconductivity may be calculated from this for-mula by setting ω = ω = − ω and the numerical resultsare shown in Figure 8. Again, we used unit cells ineach lattice direction and Chebyshev moments andcompare the results with the ones obtained by integrat-ing in k -space, just like in the previous subsection. Forconvenience, we define the constant σ = e a/ t (cid:126) [10] interms of the hopping integral t and the lattice parameter a . This particular example benefits considerably from thecancellation of the most complicated objects that neededto be calculated. In appendix VIII A we present an ex-ample with less symmetry that confirms the completeagreement between our method and the k -space formal-ism. C. Photoconductivity of h-BN with Andersondisorder
Our formalism does not rely on translation invariance,and so may be used to study disordered systems. Toshow this, we now introduce to h-BN a simple model fordisorder by letting each atomic site have a random localenergy taken from a uniform distribution [ − W/ , W/ (Anderson disorder [1]): H W = (cid:88) R (cid:88) σ W σ ( R ) c † σ ( R ) c σ ( R ) where R is the position of the unit cell and σ labels the Figure 8: Second-order yyy photoconductivity for h-BN. Hop-ping parameter: t = 2 . eV, temperature: T = 0 K, chemicalpotential: µ = 0 eV, gap ∆ = 7 . eV broadening parameter: λ = 39 meV, number of Chebyshev moments used: M = 2048 ,lattice size: L = 4096 × . The imaginary part disappearsafter the photoconductivity is properly symmetrized.Figure 9: Photoconductivity of hexagonal Boron Nitride inthe presence of Anderson disorder of varying strength W anda broadening parameter of λ = 23 meV. The parameters arethe same as for Figure 8 except for the number of polynomials,which is M = 512 . The dashed lines represent the imaginarypart of the conductivity. atoms inside each unit cell. The presence of disorder isexpected to smooth out the sharp features of the pho-toconductivity. As disorder increases, we should see adecrease in conductivity due to Anderson localization.This is the exact behavior that is seen in Figure 9 wherewe plot the photoconductivity of h-BN in the presence ofAnderson disorder of varying strength. Some fluctuationsexist at the features, which are expected to disappear asthe system size increases.As expected, the introduction of Anderson disorder produces a broadening of the sharp features of the non-linear optical conductivity. This broadening also meansthat there will be a larger response to the external electricfield at frequencies smaller than the gap.The large oscillations near the origin reveal somethinginteresting about the numerical details of our formalism.Eq. 27 is comprised of a complicated sum of severalterms. Individually, some of these terms may be verylarge, but there may be cancellations among them. Foreach of these terms, the Chebyshev expansion is exact inthe limit of infinite polynomials. For a finite number ofpolynomials, there will be slight differences between theexact result and the expansion, and if the exact resultis very large, this difference will be considerable. It ishighly unlikely that this difference will be the same foreach term, and so their sum may not cancel out in theend. This is the typical behavior at the lower-frequencyregime that is related to the singularities that plaguedthe velocity gauge approach. This has been discussedsince the early work of Sipe and challenged for a longtime the equivalence between the velocity and lengthgauges [17, 21, 22]. This effect could be fully mitigated bygreatly increasing the number of polynomials, but herewe are interested in the finite frequency behavior. D. Second-harmonic generation of h-BN withvacancies
In realistic samples, vacancies and impurities may ex-ist due to imperfections in the fabrication process, aswell as other more complex structural defects. In thissection, we show that our method allows us to obtainthe second-harmonic generation (SHG) of a system withstructural disorder. Using eq. 27, we show in Figure10 the effect of vacancies of varying concentration in theSHG of h-BN. Unlike Anderson disorder, the addition ofvacancies to the system does not change the gap. Theirmost noticeable effect is to flatten the features of thesecond-harmonic generation. As discussed in the previ-ous section, the lower frequencies are dominated by os-cillations and would require many more polynomials tofully converge. Therefore, we omit that region and onlyrepresent the remaining regions, which have already con-verged within the desired accuracy.
E. Considerations on convergence and accuracy
In this section, we briefly discuss some convergenceproperties of our method. For a more thorough discus-sion, see [23]. The convergence to the exact value de-pends on several factors:1. Spectral methods rely on the self-averaging proper-ties of random vectors, yielding an associated vari-ance. The error bar decreases as √ N R N , where N R Figure 10: Second-harmonic generation in h-BN for a varyingconcentration of vacancies and λ = 2 . meV. The blue (red)curves represent the real (imaginary) part of the conductivity.The darker curves have a higher concentration of vacancies.System size: L = 2048 , number of polynomials: M = 512 .All the other parameters are the same as in Figure 9. is the number of random vectors and N is the sizeof the sample.2. In the thermodynamic limit of an infinite lattice,the spectrum becomes continuous and so we expectthe conductivity curve to be smooth. However, thesystems used in simulations are finite and so havea typical energy level spacing, which we denote by δε . This has important consequences for the res-olution. Details characterized by a smaller energyscale than that of δε are meaningless because theycannot be distinguished from the contribution ofindividual energy levels. The maximum resolutionis therefore limited by the energy level spacing. Forour concrete examples with the honeycomb lattice,we use δε = 3 πt/L , the energy level spacing at theDirac point in Graphene for a system of linear di-mension L .3. The resolution may be controlled through λ , thebroadening parameter of the Green’s functions.Energy differences smaller than λ become indistin-guishable from one another. On the one hand, asmall λ is required in order to resolve the sharp fea-tures of the curve accurately. On the other hand,when λ (cid:46) δε , the discrete nature of the spectrumstarts to become visible through the roughness ofthe curve. For sufficiently small λ , the expectedsharp features of the curve become indistinguish-able from the contributions of the individual energylevels. If these issues are not solved, they becomea major source of systematic error in the final re-sults. Therefore, if we want to see the expected thermodynamic limit, we have to ensure λ (cid:38) δε .In Figure 11, the yy optical conductivity of Graphene isrepresented for several values of λ . In this example, δε =5 . meV. As λ is decreased, the curve becomes sharper,but when λ = 2 . meV the discreteness of the spectrumstarts to become noticeable through the roughness of thecurve. It is starting to diverge from the expected smoothcurve of the thermodynamic limit.In the lower inset, we study the convergence as a func-tion of the number of polynomials at (cid:126) ω = 4 . eV, aregion of rapidly changing conductivity. The smaller the λ , the more polynomials are required in order to have afully converged result. Within the accuracy δσ/σ (cid:39) . ,all the curves have already converged at . × polyno-mials. These calculations were repeated for several dif-ferent initial random vectors. In the plot we show onlyone of these calculations. The error bar associated withthe random vectors is too small to be distinguished fromthe curves themselves.In the upper inset, we do the same thing, but now in avery small region around (cid:126) ω = 2 . eV, a region of slowlyincreasing conductivity. The plot shows three sets ofcurves with different colors. Inside each set, we representa collection of frequencies, ranging from (cid:126) ω = 2 . eVto (cid:126) ω = 2 . eV. The darker curves correspond tohigher frequencies. The main graph shows that all thesecurves have converged to the same value in a region ofslowly increasing conductivity. The inset, however, showsa different picture. The red ( λ = 23 meV) and black( λ = 230 meV) sets of curves show a variation consistentwith the expected increasing conductivity. If one zoomsin to those sets of curves, it is possible to check thatthey are indeed increasing in value as ω increases. Thegreen curve ( λ = 2 . meV) is not only changing in a scalemuch larger than expected, but it is also decreasing. Thisvariation comes from the individual contribution of theenergy levels, not from features of the conductivity andis therefore artificial. Within the accuracy δσ/σ (cid:39) − each of these curves has completely converged at . × polynomials but this level of accuracy is meaningless for λ = 2 . meV. The error bars are not shown for clar-ity, but their values are the following: at λ = 230 meV, δσ/σ = 10 − ; at λ = 23 meV, δσ/σ = 3 × − ; at λ = 2 . meV, δσ/σ = 5 × − . At this scale, the errorbars are comparable to the variation due to the numberof polynomials and to the value of λ .These frequencies were chosen to compare the conduc-tivity in a place where it is expected to converge quicklyand another where it is expected to converge slowly.Looking at these graphs, it is possible to estimate howmany polynomials are required to converge to the finalvalue of the conductivity for the specified parameter λ within a given accuracy. A rough estimate of the scalingis given by N ∼ λ − .A similar analysis may be done for the second-order1 Figure 11: First-order optical yy conductivity of Graphene for M = 16384 and L = 2048 now as a function of the broaden-ing parameter λ . The remaining parameters remain the sameas for Figure 7. The solid (dashed) curves represent the real(imaginary) part of the conductivity. The legend shows thevalues for the broadening parameter. The lower inset showsthe evolution of the value of the conductivity for each λ asthe number of polynomials is increased for (cid:126) ω = 2 . eV. Theupper inset shows the same thing but for several very close fre-quencies around (cid:126) ω = 4 . eV. The darker curves correspondto higher frequencies. conductivity. We will not present it here for two reasons.Firstly, the main points of the previous paragraphs re-main the same. Secondly, we cannot do such an analysisbecause the computational cost would be tremendouslyhigher. F. Considerations on efficiency
Our formalism provides a very general framework withwhich to compute the nonlinear optical response up toany order. Once the formulas were obtained, we chose touse spectral methods to perform the computation. Thisis not always the most efficient approach: for systemswith translation invariance and periodic boundary con-ditions, we can specify the formulas for k -space and thenperform the explicit integration. Then, for a given set ofparameters (temperature, broadening, Fermi energy) thecomputation time will scale as L D N ω where L D is thenumber of points in the Brillouin zone (which is also thenumber of lattice sites), D is the dimensionality and N ω is the number of frequencies we want to compute. Foreach k and each frequency, this method comes down todiagonalizing the Bloch Hamiltonian H k , and then sum-ming over the whole set of k points. This method isextremely efficient at computing the optical conductivityat any order using the velocity gauge.Using spectral methods, the computation is split intothe calculation of the Chebyshev moments and the final matrix product of the Γ matrices with the Λ matrices.The first part is the most demanding and is indepen-dent of the parameters mentioned above. Its computa-tion time scales as L D N n +1 , where n is the order of theconductivity and N is the number of Chebyshev polyno-mials. More concretely, if we want to calculate the con-ductivity for a certain λ , using N ∼ λ − ∼ N ω , we findthat the k -space calculation scales much more favorably.If the system has no translation invariance, k -spaceintegration is no longer useful and we would need to nu-merically diagonalize the full Hamiltonian. This methodscales as L D N ω which is highly unfavorable and becauseof that we would be limited to very small systems. In thiscontext, spectral methods become the preferred choice.For the examples used in this paper, the computationof the second-order conductivity with the k -space formal-ism in a system with L = 2048 took around 2 minuteson a Xeon E5-2650 with 16 threads. In comparison, thesame computation took hours for translation-invarianth-BN with polynomials, and hours for h-BNwith Anderson disorder/vacancies and polynomials.Despite the discrepancy in computational efficiency, weknow of no other more efficient way to compute the non-linear optical conductivity for disordered systems. VI. CONCLUSIONS
We developed an out-of-equilibrium expansion of thenon-linear optical response of non-interacting systemsusing the Keldysh formalism and expressed everythingin terms of traces of operators. This provides a basis-independent expression for the linear and non-linear op-tical conductivities. This drops the requirement of trans-lation invariance and allows us to include magnetic fieldsand disorder in our tight-binding calculations. We alsoprovide a diagrammatic representation of this expansion,which makes it a very straightforward process to obtainthose expressions.The expressions for the non-linear conductivities arecalculated numerically with resort to an expansion inChebyshev polynomials and a stochastic evaluation ofthe trace, in close resemblance to the Kernel Polyno-mial Method (KPM). We provide the mapping that takesthe aforementioned expressions and converts them to anumerically-suited object to be calculated with spectralmethods. This is only possible because of the careful wayin which these expressions were constructed in the firstplace.We built an open-source software that is able to calcu-late the first- and second-order optical conductivities ofvery large 2D tight-binding systems ( atoms) withdisorder and magnetic fields. This software is usedto obtain the first- and second-order conductivities ofGraphene and hexagonal Boron Nitride. These samequantities were calculated independently with the usual2integration in k -space and the results are in great agree-ment, proving the validity of our method. Finally, weshow two examples with disorder that cannot be treatedunder the usual k -space formalism: h-BN with Andersondisorder and vacancies.We briefly discuss the convergence properties of thismethod by analyzing how the curves change as the res-olution is increased. The resolution is controlled by λ ,the broadening parameter of the Green’s functions andis limited by δε , the energy level spacing of the system.The expected thermodynamic limit is obtained by de-creasing δε (increasing the system size) while ensuring δε (cid:46) λ .For systems with translation invariance, the k -spaceintegration is very quick and is preferred over ourmethod. If the systems do not have this property, thismethod becomes the most efficient one to calculate thesecond-order conductivity with disorder. This paperserves as a proof of concept and the effects of realisticdisorder on the nonlinear optical properties of 2D mate-rials will be explored in a future paper. VII. ACKNOWLEDGEMENTS
The work of SMJ is supported by Fundaçãopara a Ciência e Tecnologia (FCT) under the grantPD/BD/142798/2018. The authors acknowledge financ-ing of Fundação da Ciência e Tecnologia, of COM-PETE 2020 program in FEDER component (EuropeanUnion), through projects POCI-01-0145-FEDER-02888and UID/FIS/04650/2013. The authors also acknowl-edge financial support from Fundação para a Ciên-cia e Tecnologia, Portugal, through national funds,co-financed by COMPETE-FEDER (grant M-ERA-NET2/0002/2016 – UltraGraf) under the PartnershipAgreement PT2020. We would also like to thank JMBLopes dos Santos, GB Ventura, DJ Passos, B Amorimand D Parker for their helpful comments and suggestionsregarding this paper.
VIII. APPENDIX AA. Sublattice displacement
The calculation of the photoconductivity of h-BNwas very efficient due to the cancellation of the three-dimensional Γ matrices. In this appendix, we provide anextra example, which does not benefit from that property.By changing the relative position of the two sublatticesin h-BN, we are able to obtain non-zero values in all the Γ matrices, which enables us to test the remainder ofthe formula. All the hopping parameters in this systemare exactly the same as in regular h-BN. The only differ-ence is in the distance between atoms, which changes the a a 𝛿 𝛿 𝛿 B A xy Figure 12: Displaced honeycomb lattice and choice of primi-tive vectors. velocity operators while keeping the Hamiltonian intact(See Figure 12).The primitive vectors are identical, but the nearest-neighbor vectors are different: δ = a (cid:32) √ , − (cid:33) δ = a (cid:18) , (cid:19) δ = a (cid:32) − √ , − (cid:33) . One of the sublattices was translated in the y direc-tion by a/ . The second-order xxx conductivity remainszero, but now the xxy photoconductivity is no longerzero and can be seen in Figure 13. The lattice size andnumber of polynomials used was reduced to and respectively, due to the greatly increased computationalcost. At lower frequencies, the results start to divergebecause there are not enough polynomials to resolve thisregion. The results are in great agreement with the onesobtained by k -space integration. The small oscillationsin the imaginary part are expected to disappear and thenumber of polynomials is increased. [1] P. W. Anderson. 1958. Absence of Diffusion in CertainRandom Lattices. Phys. Rev.
109 (Mar 1958), 1492–1505.Issue 5. https://doi.org/10.1103/PhysRev.109.1492 [2] M Bass, PA Franken, JF Ward, and G Weinreich. 1962.Optical rectification.
Physical Review Letters
9, 11(1962), 446.[3] Francesco Bonaccorso, Zicai Sun, Tawfique Hasan, andAC Ferrari. 2010. Graphene Photonics and Optoelec-tronics. 4 (06 2010).[4] Richard Henry Dalitz, Rudolf Ernst Peierls, and J Peifer.1997.
Selected scientific papers of Sir Rudolf Peierls: withcommentary . World Scientific, Singapore. https://cds.cern.ch/record/317707 Figure 13: Second-order xxy conductivity for displaced h-BN. Hopping parameter: t = 2 . eV, temperature: T = 0 K,chemical potential: µ = 0 eV, gap ∆ = 7 . eV broadeningparameter: λ = 39 meV, number of Chebyshev moments used: M = 512 , lattice size: L = 1024 × .[5] Aires Ferreira and Eduardo R. Mucciolo. 2015. CriticalDelocalization of Chiral Zero Energy Modes in Graphene. Phys. Rev. Lett.
115 (Aug 2015), 106601. Issue 10. https://doi.org/10.1103/PhysRevLett.115.106601 [6] RP Feynman. 1963. RP Feynman and FL Vernon, Ann.Phys.(NY) 24, 118 (1963).
Ann. Phys.(NY)
24 (1963),118.[7] P. A. Franken, A. E. Hill, C. W. Peters, and G. Wein-reich. 1961. Generation of Optical Harmonics.
Phys.Rev. Lett. https://doi.org/10.1103/PhysRevLett.7.118 [8] Walter Gautschi. 1970. On the construction of Gaussianquadrature rules from modified moments.
Math. Comp.
24, 110 (1970), 245–260.[9] M.M. Glazov and S.D. Ganichev. 2014. High frequencyelectric field induced nonlinear effects in graphene.
Physics Reports https://doi.org/10.1016/j.physrep.2013.10.003
High frequencyelectric field induced nonlinear effects in graphene.[10] F Hipolito, Thomas G Pedersen, and Vitor M Pereira.2016. Nonlinear photocurrents in two-dimensional sys-tems based on graphene and boron nitride.
Physical Re-view B
94, 4 (2016), 045434.[11] Radi A Jishi. 2013.
Feynman diagram techniques in con- densed matter physics . Cambridge University Press.[12] Leonid Veniaminovich Keldysh. 1964. Diagram techniquefor nonequilibrium processes.
Zh. Eksp. Teor. Fiz.
Journal ofPhysics: Condensed Matter
20, 38 (2008), 384204. http://stacks.iop.org/0953-8984/20/i=38/a=384204 [15] GHC New and JF Ward. 1967. Optical third-harmonicgeneration in gases.
Physical Review Letters
19, 10(1967), 556.[16] Daniel E Parker, Takahiro Morimoto, Joseph Orenstein,and Joel E Moore. 2018. Diagrammatic approach to non-linear optical response with application to Weyl semimet-als. arXiv preprint arXiv:1807.09285 (2018).[17] DJ Passos, GB Ventura, JM Viana Parente Lopes,JMB Lopes dos Santos, and NMR Peres. 2018. Non-linear optical responses of crystalline systems: Resultsfrom a velocity gauge analysis.
Physical Review B
97, 23(2018), 235446.[18] Rudolph Peierls. 1933. Zur Theorie des Diamagnetismusvon Leitungselektronen.
Zeitschrift für Physik
80, 11-12(1933), 763–791.[19] RA Sack and AF Donovan. 1971. An algorithm for Gaus-sian quadrature given modified moments.
Numer. Math.
18, 5 (1971), 465–478.[20] Julian Schwinger. 1961. Brownian motion of a quantumoscillator.
J. Math. Phys.
2, 3 (1961), 407–432.[21] J. E. Sipe and Ed Ghahramani. 1993. Nonlinear opticalresponse of semiconductors in the independent-particleapproximation.
Phys. Rev. B
48 (Oct 1993), 11705–11722. Issue 16. https://doi.org/10.1103/PhysRevB.48.11705 [22] G. B. Ventura, D. J. Passos, J. M. B. Lopes dos Santos,J. M. Viana Parente Lopes, and N. M. R. Peres. 2017.Gauge covariances and nonlinear optical responses.
Phys.Rev. B
96 (Jul 2017), 035431. Issue 3. https://doi.org/10.1103/PhysRevB.96.035431 [23] Alexander Weisse, Gerhard Wellein, Andreas Alvermann,and Holger Fehske. 2006. The kernel polynomial method.
Reviews of modern physics
78, 1 (2006), 275.[24] Fourier convention: f ( t ) = (2 π ) − ´ dωe − iωt ˜ f ( ω ))