{'fraction_non_alphanumeric': 0.038500530472710126, 'fraction_numerical': 0.008770482140752092, 'mean_word_length': 4.383424292422896, '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': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Study of Simulation Method of Time Evolution of Atomic and Molecular Systems by Quantum Electrodynamics

We discuss a method to follow step-by-step time evolution of atomic and molecular systems based on QED (Quantum Electrodynamics). Our strategy includes expanding the electron field operator by localized wavepackets to define creation and annihilation operators and following the time evolution using the equations of motion of the field operator in the Heisenberg picture. We first derive a time evolution equation for the excitation operator, the product of two creation or annihilation operators, which is necessary for constructing operators of physical quantities such as the electronic charge density operator. We then describe our approximation methods to obtain time differential equations of the electronic density matrix, which is defined as the expectation value of the excitation operator. By solving the equations numerically, we show"electron-positron oscillations", the fluctuations originated from virtual electron-positron pair creations and annihilations, appear in the charge density of a hydrogen atom and molecule. We also show that the period of the electron-positron oscillations becomes shorter by including the self-energy process, in which the electron emits a photon and then absorbs it again, and it can be interpreted as the increase in the electron mass due to the self-energy.


I. INTRODUCTION
The elementary processes of almost every phenomenon in condensed matter physics and chemistry can be regarded as time evolution of a system which consists of interacting charged particles and photons. The physical theories which describe such a system are electromagnetism and quantum mechanics, and their unified theory has already been constructed before the middle of last century as the quantum electrodynamics (QED), taking the form of quantum field theory. The QED is the most stringently tested theory of physics and its precision is confirmed by many types of experiments of elementary particle physics. It can be considered as the most successful fundamental physical theory we have.
However, there has been only limited use of QED in the fields such as condensed matter physics and quantum chemistry. There are a lot of works to simulate the time evolution of the quantum systems involving light and matter by using the time-dependent Schrödinger, Dirac, or DFT (density functional theory) equation with classical electromagnetic fields (i.e. semi-classical approximation) [1][2][3][4], but both matter and light are not treated as quantized fields. In some simulations, the quantum photon field is used but the matter part and its interaction with the photon are much more simplified than QED [5][6][7], or the quantum Dirac field is used for electrons but the interaction is not given by the photon field [8,9]. In atomic physics and quantum chemistry, QED is only regarded as a small correction to the Coulomb potential which appears in the time-independent Dirac equation [10][11][12], and not considered in dynamical situations.
It is true that the approximations above are valid for a broad range of systems so far, but the atomic, molecular, and optical physics experiments are in rapid progress recently.
It is now possible to measure, fabricate, and control structures on the atomic and molecular scale with the advances in nanotechnology. A single photon and a single electron spin can be measured and manipulated by the current technology of photonics and spintronics [13][14][15][16].
Moreover, developments in the laser science have made it possible to observe a phenomenon at the time scale of femtosecond to attosecond time scales [17,18]. Considering such progress in the experiments of smaller space-time scales and the most fundamental particle properties, it is important for the theoretical side to develop a simulation method based on as an elementary theory as possible. This is the reason why we try to formulate a time evolution simulation method for atomic and molecular systems, which consist of electrons, atomic nuclei and photons, closely following QED in the form of quantum field theory.
To achieve our goal to develop a method to simulate time evolution of atomic and molecular systems by QED, we have to overcome some issues which do not appear in the ordinary QED. When we say the ordinary QED, we mean that it is a relativistic (Lorentz invariant) quantum field theory and the calculation method of the scattering amplitude is performed by the covariant perturbation theory [19,20]. It implies that the only transition between infinite past ("in-state") and infinite future ("out-state") is concerned and the quantum fields in those states ("asymptotic states") can be treated like non-interacting fields. This in turn makes it possible to apply the perturbative approach by taking the asymptotic states (where fields are described by non-interacting Hamiltonian) as the unperturbed states and interaction as the perturbation. Also, the perturbative method is established in a quite systematic way, owing to the Lorentz invariance of the theory. Although such a method yields physical quantities which can be precisely compared with experiments of particle physics in particular, it is not sufficient for our interests. For simulating time evolution of atomic and molecular systems, there are mainly three issues we have to face which are not simultaneously concerned in the ordinary QED: (i) there are atomic nuclei, which are non-relativistic and not elementary particles, (ii) the matter particles in those systems are in bound states, (iii) we would like to follow the finite time evolution of the systems step by step.
In fact, there have already been methods which can partly treat these three points within the framework of quantum field theory. First, as for (i), there are effective field theories in which nucleons are treated as quantized fields [21,22], but the computation is within the scattering theory, so they are not suited to our needs of (ii) and (iii). In principle, since the nucleons consist of quarks and gluons which are described by quantum chromodynamics (QCD), they can be treated in a quantum field theoretic way, and the method of lattice field theory has been well developed to perform non-perturbative calculation of QCD, known as lattice QCD calculation. However, it takes too much computational time even to describe a nucleon as a bound state of quarks, and adding QED to study atoms is much more unpractical. Other shortcomings include that the current lattice field theory is developed only to describes an equilibrium state and time evolution cannot be treated. Also, in contrast to lattice QCD, lattice QED has a problem that, being U(1) gauge theory which does not exhibit asymptotic freedom, the continuum limit of lattice spacing cannot be taken. Second, as for (ii), a well-known technique is the Bethe-Salpeter equation [23]. In Ref. [24], the Bethe-Salpeter equation and various other techniques to describe bound states in QED are reviewed, but those incorporating (i) or (iii) are not found. Finally, as for (iii), some formalisms are known to describe a non-equilibrium state by quantum field theory, such as closed time path (CTP) formalism [25,26] and thermo field dynamics [27]. In particular, CTP formalism is applied to gauge theory including QED and QCD, but currently not for bound state problems. This is because the CTP formalism uses systematic perturbative expansion based on the interaction picture, but we cannot divide the QED hamiltonian into unperturbed and interaction parts when we consider bound states. Therefore, to simulate time evolution of atomic and molecular systems based on QED, we cannot just use preexisting formalisms as mentioned above. We briefly explain our proposals [28,29] to cope with these issues in the following. First of all, as for (i), we add the atomic nuclei degree of freedom as Schrödinger fields, and the interaction with the photon field is determined from U(1) gauge symmetry as usual [28]. Although this deprives the theory of the Lorentz invariance, since we try to follow the finite time evolution of the system, this would not be a crucial problem. Next, as for (ii), in order to describe the bound state, we expand the matter field operators by localized wavepackets, not by the usual plane waves, and define the creation and annihilation operators [29]. In the case of the Dirac field operator for electrons, for example, we may adopt the stationary solutions of the Dirac equation under the existence of external electrostatic field, as the expansion functions. This is similar to the Furry picture [30,31], but, as described below, we do not assume the time-dependence of the operators which is determined by the energy eigenvalues. Finally, as for (iii), we follow the time evolution using the equations of motion of the field operator in the Heisenberg picture [29]. As is mentioned earlier, since we do not have well-defined division between unperturbed and interaction parts of the QED hamiltonian for the bound state, we cannot work in the interaction picture.
The setups described above determine the evolution equations of field operators, and those of creation and annihilation operators, but we have to make further approximations and assumptions to obtain the evolution equations for the expectation value of physical quantities. In Ref. [29], we have studied the time evolution of one of the most basic physical quantity operators, the electronic charge density operator, and have discussed approximation methods to obtain the time evolution of its expectation value. The charge density operator is expressed by the product of two creation or annihilation operators, which is called an excitation operator, but since the time derivative of creation and annihilation operators contains more than one of these operators, the time differential equation of the excitation operator is not closed. In other words, differentiating the excitation operator yields operators which cannot be expressed by the excitation operator. Such a problem is generic for interacting quantum field theories [25], and not special to our approach. In Ref. [29], we have introduced several approximations to obtain the time evolution equation for the expectation value of the excitation operator, which is called a density matrix, and numerically solved the time evolution of the density matrix. The time evolution of the expectation value of the charge density operator has been obtained by multiplying the density matrix by the expansion functions of the field operator. Then, we have found that the time evolution of the charge density of a hydrogen atom exhibits very rapid oscillations of the period ≈ 1.7×10 −4 a.u. In the present paper, we improve one of the approximations employed in Ref. [29] with respect to the terms which include the photon creation and annihilation operators. In our formalism, the time derivative of the electron excitation operator has the terms which consist of two creation or annihilation operators sandwiching a photon creation or annihilation operator (we call this type of operator by "ê †âê -type operator" for short). In Ref. [29], when we take the expectation value of these terms, we have factorized the terms into the expectation value of the excitation operator and that of the photon creation or annihilation operator. After this approximation, these terms give finite contribution only when an initial photon state is a coherent state, which is an eigenstate of the photon annihilation operator.
In particular, they vanish for the photon vacuum state. To go one step further, we, in this paper, do not perform the above factorization, and solve the time evolution equation of thê e †âê -type operators simultaneously with that of the excitation operator. As we show in a later section, this procedure corresponds to counting the self-energy process of the electron (the electron emits a photon and then absorbs it again), and, consistently, it gives non-zero contribution even when the initial state is the photon vacuum state. This paper is organized as follows. In Sec. II, we derive the time evolution equations of the quantum operators. After we describe how we expand field operators and define creation and annihilation operators, we show time derivative of the creation and annihilation operators from the equations of motion of the quantum fields. Then, we derive the time evolution equation of the excitation operator andê †âê -type operator. In Sec. III, we derive time evolution equations of the density matrix by taking the expectation value of the time derivative of the excitation operator. We explain our approximation methods to derive closed sets of time evolution equations. In Sec. IV, we show the results of numerical computation of these equations for a hydrogen atom and molecule. Finally, Sec. V is devoted to our conclusion. The notations and conventions follow those in Refs. [28,29,32,33].

II. EVOLUTION EQUATIONS FOR QUANTUM OPERATORS
We would like to describe atoms and molecules as systems which consist of electrons, atomic nuclei and photons. In this paper, we work in the Born-Oppenheimer (BO) approximation, in which positions of the atomic nuclei are fixed. Then, only electrons and photons are described as quantum field operators, and atomic nuclei contribute to the charge density as delta functions. The source of the photon field can be both electrons and atomic nuclei.
The quantum field operators which appear in this paper are the four-component Dirac field operatorψ(x) for the electron, and the U(1) gauge fieldÂ µ (x) for the photon. (Incidentally, the approach to treat the electron as the two-component Schrödinger field [34] is also developed in our group. See Refs. [35][36][37] for details.) Although there are some overlaps in this section with the contents in Ref. [29], we reproduce them in the reorganized form for the convenience of the readers.

A. Definitions of creation and annihilation operators
Our expansion of the Dirac field operator iŝ where ψ n + ( r) and ψ n − ( r) are respectively the n-th electron and positron solutions of the four-component Dirac-Hartree-Fock equation under the existence of external electrostatic field, and they form an orthonormal basis set as d 3 r ψ † n a ( r)ψ m b ( r) = δ nm δ ab . In Eq. (1), N D is the number of the electron expansion functions, which is same as one of the positron expansion functions. Considering the Kramers pair, N D usually equals to twice the number of basis functions used to solve the Dirac equation. We note that N D has to be infinite to make the expansion function set complete, but this is not available in numerical calculation.
In practice, we use a finite set which only spans a certain region of the complete space. We should interpret thus obtained results as phenomena within the subspace spanned by the finite set of expansion functions. We may use any localized wavepackets for the expansion functions, but the above choice is convenient because we can obtain such functions easily by the publicly available code like DIRAC [38]. The creation and annihilation operators are defined as the coefficients of the expansion functions and carry the time dependence.
In our notation,ê n + is the electron annihilation operator andê n − is the positron creation operator. Accordingly,ê † n + is the electron creation operator andê † n − is the positron annihilation operator. In the literature, a creation operator is usually expressed by an operator with a dagger as superscript, but, note that, in our notation, the positron creation operator does not carry a dagger whereas the positron annihilation operator does. We adopt this notation to make compact the expression of the expansion of Eq. (1) and equations below.
The equal-time anti-commutation relation can be written as ê n a (t),ê † m b (t) = δ nm δ ab , and anti-commutators of other combinations are zero.
For later convenience, we here define the electronic excitation operator, which is formed from two creation or annihilation operators aŝ Then, the electronic charge density operator and current density operators are written by the excitation operator asρ where we define and Z e = −1. The total charge density operator have contribution from both electrons and atomic nuclei,ρ(x) =ρ e (x) + Nn a=1ρ a (x), whereρ a (x) is the atomic nuclear charge density operator, a denotes the type of atomic nucleus and we assume N n types of atomic nuclei where Z a is the nucleus a's atomic number, and R a should be understood as the direct sum of the position of each nucleus of type a (see Appendix of Ref. [29] for details). As for the total charge current density operator, since there is no contribution from the atomic nuclei under the BO approximation, we haveˆ j(x) =ˆ j e (x). Although we do not compute the energy of the system in this paper, it may be instructive to show the QED Hamiltonian operator. This is shown in Appendix A.
As for the photon field operator, we use the integral equation form of the equation of motion. Namely, we use the fact that the formal solutions of the inhomogeneous Maxwell equations are known from the classical electrodynamics [5,39,40], and expressÂ µ (x) by such solutions [28,29]. Similar technique is used to derive the so-called Yang-Feldman equation, which is originally introduced to discuss the S-matrix of QED in Heisenberg picture [41].
We adopt the Coulomb gauge, ∇ ·ˆ A(x) = 0. Then, the scalar potential is given bŷ and the vector potential byˆ where the first term is the quantized free radiation field in the Coulomb gauge [5,42] and the second term is the retarded potential.ˆ A rad (ct, r) is expressed aŝ whereâ pσ (â † pσ ) are the annihilation (creation) operator of the photon with momentum p and helicity σ, and e is the polarization vector. They satisfy the commutation relation [â pσ ,â † qτ ] = δ( p− q)δ στ , and commutators of the other combinations are zero. Our convention for e( p, σ) is found in Ref. [29].ˆ A A (ct, r) is expressed aŝ where we define the retarded time u = t − | r − s|/c, andˆ j T (x) is the transversal component of the current. Writing explicitly, where As we mentioned in Sec. I, we do not consider the usual asymptotic states of the field operator, in which the state becomes non-interacting as t → −∞ asymptotically. We assume that the interaction is absent for t < 0 and the free part expressed above byÂ rad,µ (x) is considered to realize at t = 0. This, in particular, implies the initial condition and enables us to rewrite the retarded potential as [34] A A (ct, r) = 1 To derive this, we use the following formulae for the delta function δ( . This form may be convenient for the numerical calculation since the retarded time is eliminated, but at the cost of increasing the dimension of integration. In passing, some comments on our assumption at the initial time (t = 0) may be in order. Following Ref. [43], we assume that the parameters and fields which appear in the equations of motion have been renormalized at t = 0. Specifically, the renormalization has been performed in a standard manner as where "B" in the subscripts denotes a bare parameter or a bare field, and δm e , Z 2 and Z 3 are the renormalization constants. However, it does not mean that we do not need renormalization at t > 0. In fact, as is discussed later in this paper, the electron mass is shown to be increased by including the self-energy process in our time evolution simulation.

B. Time evolution of annihilation and creation operators
The time derivative ofê n a is given by substituting the expansion Eq. (1) into the Dirac field equation, multiplying by ψ † n a ( r), integrating over r, and using the orthonormality condition [29]. This leads to where the various coefficient matrices are defined as follows. In the first term, we define where T n a m b is the electronic kinetic energy integral, M n a m b is the electronic mass energy integral, and V n a m b ( R) is the nuclear attraction integral, respectively defined as They originate respectively in the kinetic term, mass term, and the nuclear charge contribution to the scalar potential term.
In the second term of Eq. (15), we define the electronic repulsion integral This term originates in the electronic charge contribution to the scalar potential term. The third and forth terms of Eq. (15) originate in the vector potential term. The former comes from the retarded partˆ A A (10) and the latter from the free radiation partˆ A rad (9). Finally, in the forth term of Eq. (15), we define C. Time evolution of excitation andê †âê -type operators In this paper, we are interested in the time evolution of electronic charge density, and its operatorρ e (x) is expressed by the excitation operator as Eq. (3). Since the excitation operator carries every time dependence ofρ e (x), what we need to know is the time evolution equation of the excitation operator.
The time derivative of the excitation operatorÊ n a m b (2) can be written as where we defineÔ Using Eq. (15), this can be readily obtained as where we have defined to make the expression shorter. For later convenience, we call the second term of Eq. (24) "four-electron term", the third term "retarded potential term", and the forth term "radiation term".
Next, we consider the time derivative of theê †âê -type operators. Since the time derivative pσê n a is known by taking the Hermite conjugate of that ofê † n aâ pσêm b , we only show the latter below. This can be expressed as where we defineP Using Eq. (15), we obtain and We note that the last term ofP n a pσm b includes the operator product of the formê † n aâ pσâ † qτê r e . This combination of the operators could express a self-energy process, in which the electron emits a photon and then absorbs it again. The effect of the self-energy process will be discussed and numerically demonstrated in Sec. IV C.

A. Definition of density matrix
We begin by introducing notations regarding expectation values. We denote the expectation value of the excitation operator with respect to the Heisenberg initial ket |Φ by and call this quantity the density matrix. Below, we sometimes write just brackets around the operator to denote the expectation value with respect to |Φ , namely, · · · ≡ Φ| · · · |Φ where · · · stands for some operators. Since (Ê n a m b ) † =ê † m bên a =Ê m b n a , taking the expectation value yields (E n a m b ) * = E m b n a , showing that the density matrix is a Hermite matrix.
Similarly, we define the expectation value of Eq. (23) as O n a m b ≡ Φ|Ô n a m b |Φ . Then, since which is the time evolution equation for the density matrix.
Using the density matrix, the expectation value of the electronic charge density operator, Eq. (3), can be written as However, this quantity has a non-zero value when the Heisenberg initial ket is the vacuum, where we define E 0 n a m b ≡ 0|Ê n a m b |0 .

B. Four-electron term and retarded potential term
The four-electron term consists of four electron creation or annihilation operators. We approximate this term by This decomposition holds exactly for t = 0 and the approximation is motivated by this fact.
In our previous paper Ref. [29], we have only used the first term, and the second term, which describes the exchange effect, has been omitted.
As for the retarded potential term, sinceĵ k T is expressed as Eq. (11), it seems to have the same structure as the four-electron term. However, sinceĵ k T is computed at the retarded time and contains a time derivative of the excitation operator, the approximation like Eq. (35) is not applicable. Thus, we approximate this term by just replacingĵ k T by its expectation value as where we define J k T ≡ Φ|ĵ k T |Φ . With this approximation, the two space integrations over r and s in the expectation value of the retarded term becomes that of c-number. In this six-dimensional integration, the retarded time u depends both on r and s, so the integration has to be done numerically, which is not practical. To avoid the six-dimensional numerical integration, we rewrite the expression into the form in which the retarded time is eliminated. This is done by the same process as the one used to derive Eq. (14). Then, for the expectation value of the retarded term, we have where we define with Here, we put [E,Ė] after I j T for the notation of the integral (38) in order to emphasize that it depends on the density matrix and its time derivative at times earlier than t. Since the functions appearing in Eqs. (39) and (40) are defined by Eqs. (6) and (12), I jj and I jE are four-center integrals and analytic formulae can be obtained when the expansion functions of the Dirac field operator are gaussian functions. We show the integral formulae in the Appendix B.

C. Radiation term
In this section, we describe two approximation methods for the expectation value of the radiation term. We first describe the simplest approximation method, which is same as the one adopted in Ref. [29]. In this approximation, we use ê † n aâ pσêr e ≈ â pσ ê † n aê r e = â pσ E n a r e , and, similarly, ê † n aâ † pσê r e ≈ â pσ * E n a r e . This factorization holds exactly for t = 0, whenê n a andâ pσ commutes, and the approximation is motivated by this fact. We note that this term may give finite contribution only when the initial photon state is a coherent state, which is an eigenstate of the photon annihilation operator. Then, combined with the approximations described in the previous subsection, we obtain This gives us a closed differential equation for the density matrix.
In the second approximation method, we do not use the above factorization. We use evolution equation for the expectation value of theê †âê -type operators simultaneously with one for the density matrix. We first define Note that Φ|ê † n aâ † pσê m b |Φ = E * m b pσn a . We next define the expectation value of the operator P n a pσm b , Eq. (27), as P n a pσm b ≡ Φ|P n a pσm b |Φ . Similarly for the operatorQ n a pσm b , Eq. (28), we define Q n a pσm b ≡ Φ|Q n a pσm b |Φ . Then, the expectation value of Eq. (24) can be written and the time evolution of E n a pσm b can be expressed as which is obtained by taking the expectation value of Eq. (26).
As for P n a pσm b , the expectation value of the operator product in the second term of Eq. (29) is approximated to be In this approximation, we assume that the initial photon state is a number state and not a coherent state. Then, since â pσ = 0, the contribution of this term vanishes. The third term of Eq. (29) includes the retarded time and its expectation value would be approximated in the same way as the retarded potential term of Eq. (24). Following the procedure described below Eq. (36), we have the expression of the form of Eq. (37) with E n a r e (t) replaced by E n a pσr e (t). The forth term of Eq. (29) includes two four-operator terms, which are approximated to be ê † n aâ pσâ qτêr e ≈ â pσâ qτ ê † n aê r e , Since we do not consider a coherent state for the initial photon state as mentioned above, â pσâ qτ = 0 in the first equation and its contribution vanishes. In the second equation, since â pσâ † qτ = (n pσ + 1)δ (3) ( p − q)δ στ , where n pσ is the occupation number of the photon mode ( p, σ) in the initial state, it may give a non-zero contribution. Putting these approximations together, we have when we assume a number state for the photon initial state. The expectation value of Eq. (30) can be approximated in a similar manner to be where we have used â † pσâ qτ = n pσ δ (3) ( p − q)δ στ .

IV. RESULTS
In this section, we show the results of numerical solution of the time evolution equations which have been derived in the previous section. The computation is performed for a hydrogen atom and molecule using the QEDynamics code [44] developed in our group.
In Sec. IV A, we describe our setups for numerical calculation including the initial condition for the density matrix. The results of two approximation methods discussed in Sec. III C are respectively presented in Sec. IV B and Sec. IV C.

A. Setups for numerical calculation
To perform numerical calculation, we first need to determine an orthonormal set of ex- We next explain the initial condition for the density matrix. We choose the initial Heisenberg ket |Φ to be the ground states of the hydrogen atom and molecule, which are obtained by the above explained computation method. Namely, expressing |Φ = |Φ e ⊗ |Φ ph , where |Φ e is the electron part and |Φ ph is the photon part, we use |Φ e =ê † 1 + |0 for the hydrogen atom, and |Φ e =ê † 1 +ê † 1 + |0 for the hydrogen molecule. In general, the ground state of a N eelectron system in the Hartree-Fock method is expressed as |Φ e = Ne/2 i=1ê † i +ê † i + |0 when N e is even, and |Φ e =ê † when N e is odd. For later use, we here introduce the terminology "occupied" orbitals. If |Φ e containsê † i + , the i-th electron orbital is called "occupied" (i can be with or without bar), and we can write |Φ e = i=occupiedê † i + |0 . Then, using the anti-commutation relation, the initial condition for the density matrix is In particular, the vacuum expectation value at t = 0, which is needed to compute Eq. (34), atom and molecule, we report the electronic charge density at (x, y, z) = (0, 0, 1). To solve the differential equations, we use the Euler method with the time step 10 −9 a.u. In this paper, we omit the contribution from the retarded potential by setting the integral I j T [E,Ė], Eq. (38), to be zero. This integral, including numerical integration, has to be computed at every time step, and performing this straightforwardly takes too much computational time.
We shall study an effective approximation method in our future work, and just neglect it in the present work.

B. Effect of photon coherent state
In this section, we show the results when the first approximation method explained in Sec. III C is adopted. Namely, we solve the time differential equation (32) using Eq. (42).
As mentioned in the end of Sec. IV A, the third term of Eq. (42), expressing the retarded potential, is neglected. As for the photon initial Heisenberg ket, when there is no radiation field, |Φ ph = |0 , since â pσ = 0, the forth term of Eq. (42) is dropped. In fact, â pσ = 0 holds when |Φ ph is any photon number state. We may have â pσ = 0 when |Φ ph is a coherent state, and we study its effect in this section. Since we quantize the radiation field in the whole space, we consider the continuous-mode coherent state. Following the notation of Ref. [5], we denote it as |{α} . This is characterized by the eigenvalues of the photon annihilation operator of each mode ( p, σ) aŝ where α( p, σ) is called spectral amplitude [5]. In this paper, we use the delta-function type spectral amplitude and its center is chosen to be a mode ( p j , σ j ), as Such photon state corresponds to a classical oscillating electromagnetic field whose propagating direction is p j , direction of circularly polarization is σ j , and amplitude is proportional to α j . The period of the oscillations is determined by p 0 = | p j | as 2π/(p 0 c).
We first show the results when there is no initial radiation field. The case of the hydrogen atom is shown in the upper panel of Fig. 1 and the case of the hydrogen molecule is in the upper panel of Fig. 2. In these figures, the variation of the electronic charge density from its initial value is plotted. The common feature is the oscillations with very short period of about 1.7 × 10 −4 a.u. As has been argued in Ref. [29], since this is very close to the period which is determined from twice the mass of electron, 2π/(2m e c 2 ) = 1.67 × 10 −4 , it can be interpreted as the fluctuations originated from virtual electron-positron pair creations and annihilations. Hence, we call this phenomenon the "electron-positron oscillations" [29]. In Ref. [29], we show this using the hydrogen atom. In this paper, we show that the electronpositron oscillations occur similarly for the hydrogen molecule, and expect that we find them universally for any atomic and molecular systems.
We next show the results when the initial photon state is a coherent state. We choose its spectral amplitude to be the form expressed by Eq. (54) with p j /| p j | = (1, 0, 0) and σ j = +1.
We including the retarded potential may affect the electron-positron oscillations.

C. Effect of electron self-energy
In this section, we show the results when the second approximation method explained in Sec. III C is adopted. Namely, we solve time differential equations (32) and (45) using Eqs. (44), (49) and (50). As in Sec. IV B, the terms which originate from the retarded potential are neglected. As for the photon initial Heisenberg ket, we assume that there is no radiation field, |Φ ph = |0 , in this section. Therefore, the occupation number is zero for every photon mode ( p, σ), n pσ = 0. However, even in such a case, as is expressed by the factor (n pσ +1) in the third term of Eq. (49), every photon mode contributes to the radiation term. This is reasonable because, as is mentioned in the end of Sec. II C, this term comes from a self-energy process of the electron, in which the electron emits a virtual photon and then absorbs it again. This virtual photon could have any momentum.
Before showing our results, we explain here the numerical details regarding the discretization of the photon modes. We have to discretize the index p σ in E n a pσm b to perform numerical calculation. We adopt the spherical coordinate system (p 0 , θ, φ) to express p and use equally spaced grid points for each coordinate whose numbers are denoted by N p 0 , N θ and N φ respectively. One more parameter we need to specify is the maximum of p 0 , denoted by p 0 max . We first set N θ = 5 and N φ = 4, and compare the case with (p 0 max , N p 0 ) = (10, 10) and the case with (p 0 max , N p 0 ) = (20,20). Since the results do not change, we adopt (p 0 max , N p 0 ) = (10, 10) in the following. As for the choice of (N θ , N φ ), the result with (N θ , N φ ) = (7,8) is slightly different from the case with (N θ , N φ ) = (5, 4), but it is almost same as the result with (N θ , N φ ) = (11,10). Therefore, in summary, we adopt (p 0 max , N p 0 , N θ , N φ ) = (10, 10, 11, 10) as the photon-mode discretization parameters.
The results are shown in Fig. 3, the upper panel for the hydrogen atom and the lower panel for the hydrogen molecule. As in Sec. IV B, the variation of the electronic charge density from its initial value is plotted. In each panel, the solid red line shows the result without the self-energy process and the green dashed line shows one including the self-energy process. Note that computation for the case without the self-energy process is same as the one described in Sec. IV B. However, in Fig. 3, the results are multiplied by 10 4 for the hydrogen atom and by 10 2 for the hydrogen molecule in order to make easy the comparison with the case including the self-energy process. In the figure, we see the electron-positron oscillations still take place when we include the self-energy process for both hydrogen atom and molecule. However, the period of the oscillations is slightly shorter than the case without the self-energy process. We have mentioned earlier that this rapid oscillations are originated from virtual electron-positron pair creations and annihilations and their period is inversely proportional to the electron mass. Therefore, decrease in the period implies increase in the electron mass. This is reasonable because the electron mass should be increased by including the self-energy process.
We note that this "increase" in the electron mass is not the physical reality. The selfenergy of the electron is the interaction energy between the electron and the electromagnetic field which is originated from the electron itself [40,42]. (It exists for either classical or quantum electrodynamics.) Since this is something we cannot remove from the electron, the total energy including the self-energy is considered to give the observed electron mass. This is the idea of the (mass) renormalization. Although the method of the renormalization is well-established for the ordinary perturbative QED, since it is based on the notion of the asymptotic states, which exist in the infinite past and future, it is not straightforwardly applicable to our QED simulation in which finite time evolution is followed. We have succeeded in extracting the self-energy of the electron in our simulation. Our next task is how to renormalize the "increase" in the electron mass, and this will be studied in our future work.
In the end of this section, we shall make comments on the self-energy effect. Although the change in the period of the electron-positron oscillations can be interpreted as the electron mass shift due to the self-energy effect, one may wonder why the shift is not infinite as in the ordinary QED. One reason is that we have truncated the infinitely many hierarchy of time evolution equations of operators at the level ofê †âê -type operator. We would have a greater self-energy effect by considering the time evolution equation of higher order operators such as theê †ââ †ê -type operator, which appears in the last term of Eq. (29). Another reason is that we have only included localized wavepackets for the expansion functions in Eq. (1).
During the self-energy process, when the electron emits the virtual photon, the electron could be in an unbounded state as a virtual particle, but such a state cannot be expressed by our present expansion functions. In order to improve this point, we may add plane wave functions to the expansion functions to express continuum modes.

V. CONCLUSION
In this paper, we have discussed a method to follow the step-by-step time evolution of atomic and molecular systems based on QED. Our strategy includes expanding the electron field operator by localized wavepackets to define creation and annihilation operators and following the time evolution using the equations of motion of the field operator in the Heisenberg picture. Under the BO approximation, we have first derived a time evolution equation for the excitation operator, which is the product of two creation or annihilation operators. We need this operator to construct operators of physical quantities such as the electronic charge density operator. We have then described our approximation methods to obtain time differential equations of the electronic density matrix, which is defined as the expectation value of the excitation operator.
In particular, we have presented two approximation methods for the expectation value of the radiation term, which includes theê †âê -type operators. One is to factorize their expectation values into the expectation value of the excitation operator and that of the photon creation or annihilation operator, and has been used to study the effect of external oscillating electromagnetic field by setting the initial photon state as a coherent state. Another is to solve the time evolution equation of theê †âê -type operators simultaneously with that of the excitation operator, which enables us to include the self-energy effect of the electron.
By solving these equations numerically, we have shown the electron-positron oscillations appear in the charge density of a hydrogen atom and molecule, for the cases both with and without including the self-energy process. We have also shown that the period of the electron-positron oscillations becomes shorter by including the self-energy process, and it can be interpreted as the increase in the electron mass due to the self-energy.
Although the results obtained in this paper can be reasonably interpreted so far, there are many things to incorporate for establishing the time evolution simulation method of atomic and molecular systems based on QED. Two important points which are not included in the present work are computation of the retarded potential and renormalization of the electron mass. As for the retarded potential, its computation is likely to be achieved by using the gaussian integral formulae which have been derived in the Appendix, although we need efficient approximation and storage methods. As for the electron mass renormalization, we shall develop a different method from that of the ordinary QED, because our renormalization should also be performed step-by-step in time. Specifically, we may need a time-dependent renormalization factor. These issues will be addressed in our future works and incorporated in our computation code.
where the coefficient matrices are defined in Eqs. (16) and (20). Here, we have excluded some terms which are infinite constants.
By taking the expectation value of the normal-ordered product of Eq. (A2) with respect to the Heisenberg ket, we can obtain the energy of the system under the electrostatic limit, E QED,electrostatic . When the Heisenberg ket is assumed to be the one introduced in Sec. IV A, at t = 0, it gives the ordinary DHF energy as h n a n a + 1 2 where ⊕ denotes the occupied electron orbitals.
Appendix B: Molecular integral formulae for retarded potential term When we compute the retarded potential term as described in Sec. III B, we need two types of four-center integrals shown in Eqs. (39) and (40). To our knowledge, gaussian integral formulae which are needed to compute them are not seen in literature. In the present paper, the retarded potential term is neglected as explained in Sec. IV, but we shall derive the formulae in this section for convenience of future works. The formulae and derivation here are based on a method described in Ref. [45].
In the case of α = 0, we can show that where D = P − Q, (B13) and is a function defined in Ref. [45] and its recursion formula is also discussed there.
The differentiation to derive [NLM|θ|N ′ L ′ M ′ ] can be done in a straightforward manner.
As for θ jj , we can show where H n (x) is a Hermite polynomial of degree n. As for θ jE , we can show Here, we have definedR where T = α T (D 2 x + D 2 y + D 2 z ). For generating a table of allR N LM up to some maximum N + L + M, recursion relations discussed in Ref. [45] can be applied. In particular, we can The details of the recursion relations and efficient numerical techniques are found in Ref. [45].