text:
\section{Introduction}  \label{sec: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\"{o}dinger, Dirac, or DFT (density functional theory) equation with
classical electromagnetic fields ({\it i.e.} semi-classical approximation) \cite{Bauer2006, Iwasa2009, Mocken2008,Joachain},
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 \cite{Loudon,Meystre,Haroche},
or the quantum Dirac field is used for electrons but the interaction is not given by the photon field \cite{Krekora2005,Krekora2006}.
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 \cite{Dyall,Reiher,Lindgren},
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 \cite{Kimble1998,Claudon2010,Buller2010,Aharonovich2011}.
Moreover, developments in the laser science have made it possible to observe a phenomenon
at the time scale of femtosecond to attosecond time scales \cite{Drescher2002,Krausz2009}.
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 \cite{Peskin,WeinbergI}.
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 \cite{Weinberg1990,WeinbergII},
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 \cite{Salpeter1951}.
In Ref.~\cite{Grotch2002}, 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 \cite{Rammer, Calzetta} and thermo field dynamics \cite{Umezawa}.
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 \cite{Tachibana2003,Ichikawa2013} to cope with these issues in the following.
First of all, as for (i), we add the atomic nuclei degree of freedom as Schr\"{o}dinger fields, and
the interaction with the photon field is determined from $U(1)$ gauge symmetry as usual \cite{Tachibana2003}.
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 \cite{Ichikawa2013}.
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 \cite{Furry1951,LLRQT}, 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 \cite{Ichikawa2013}.
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.~\cite{Ichikawa2013}, 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 \cite{Rammer}, and not special to our approach.
In Ref.~\cite{Ichikawa2013}, 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 $\approx 1.7 \times 10^{-4}$\,a.u. ($4.1\times 10^{-21}$\,s),
which corresponds to the inverse of twice the electron mass.
This is interpreted as the fluctuations originated from the virtual electron-positron pair creations and annihilations,
showing the effect of QED, and we have designated the phenomenon as ``electron-positron oscillations".

In the present paper, we improve one of the approximations employed in Ref.~\cite{Ichikawa2013} 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 ``$\hat{e}^\dagger \hat{a} \hat{e}$-type operator" for short).
In Ref.~\cite{Ichikawa2013}, 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 the $\hat{e}^\dagger \hat{a} \hat{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.~\ref{sec:qnum}, 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 $\hat{e}^\dagger \hat{a} \hat{e}$-type operator.
In Sec.~\ref{sec:dm}, 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.~\ref{sec:results}, we show the results of numerical computation of these equations for a hydrogen atom and molecule.
Finally, Sec.~\ref{sec:conclusion} is devoted to our conclusion.
The notations and conventions follow those in Refs.~\cite{Tachibana2001,Tachibana2003,Tachibana2010,Ichikawa2013}.

\section{Evolution equations for quantum operators}  \label{sec:qnum}

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 $\hat{\psi}(x)$ for the electron,
and the $U(1)$ gauge field $\hat{A}_\mu(x)$ for the photon.
(Incidentally, the approach to treat the electron as the two-component Schr\"{o}dinger field \cite{Tachibana2013}
is also developed in our group. See Refs.~\cite{Senami2013a,Senami2013b,Senami2014} for details.)
Although there are some overlaps in this section with the contents in Ref.~\cite{Ichikawa2013},
we reproduce them in the reorganized form for the convenience of the readers.

\subsection{Definitions of creation and annihilation operators} \label{sec:qnum_defac}

Our expansion of the Dirac field operator is
\begin{eqnarray}
\hat{\psi}(ct, \vec{r}) = \sum_{n=1}^{N_D}  \sum_{a=\pm}  \hat{e}_{n^a}(t) \psi_{n^a}(\vec{r}), \label{eq:psi_expand2}
\end{eqnarray}
where $\psi_{n^+}(\vec{r})$ and $\psi_{n^-}(\vec{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
$\int d^3\vec{r}\, \psi_{n^a}^\dagger(\vec{r})\psi_{m^b}(\vec{r}) = \delta_{nm} \delta_{ab}$.
In Eq.~\eqref{eq:psi_expand2}, $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 \cite{DIRAC12}.
The creation and annihilation operators are defined as the coefficients of the expansion functions and
carry the time dependence.
In our notation, $\hat{e}_{n^+} $ is the electron annihilation operator and $\hat{e}_{n^-} $ is the positron creation operator.
Accordingly, $\hat{e}^\dagger_{n^+} $ is the electron creation operator and $\hat{e}^\dagger_{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.~\eqref{eq:psi_expand2} and equations below.
The equal-time anti-commutation relation can be written as
$\left\{ \hat{e}_{n^a}(t), \hat{e}_{m^b}^\dagger(t) \right\} = \delta_{nm} \delta_{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 as
\begin{eqnarray}
\hat{\cal E}_{n^a m^b} \equiv \hat{e}^\dagger_{n^a} \hat{e}_{m^b}.  \label{eq:op_excitation}
\end{eqnarray}
Then, the electronic charge density operator and current density operators are written by the excitation operator as
\begin{eqnarray}
\hat{\rho}_e(x) &=&  \sum_{n,m=1}^{N_D} \sum_{a,b=\pm} \rho_{n^a m^b}(\vec{r}) \hat{\cal E}_{n^a m^b}(t), \label{eq:op_rhoe} \\
\hat{j}_e^k(x) &=& \sum_{n,m=1}^{N_D} \sum_{a,b=\pm}  j^k_{n^a m^b}(\vec{r}) \hat{\cal E}_{n^a m^b}(t),  \label{eq:op_je}
 \end{eqnarray}
where we define
\begin{eqnarray}
\rho_{n^a m^b}(\vec{r}) &\equiv& (Z_e e) \psi_{n^a}^\dagger(\vec{r})  \psi_{m^b}(\vec{r}),  \\
j^k_{n^a m^b}(\vec{r}) &\equiv& Z_e\, e\,c \left[ \psi_{n^a}^\dagger(\vec{r}) \gamma^0 \gamma^k \psi_{m^b}(\vec{r})  \right], \label{eq:jk}
\end{eqnarray}
and $Z_e=-1$.
The total charge density operator have contribution from both electrons and atomic nuclei,
$\hat{\rho}(x) = \hat{\rho}_e(x) + \sum_{a=1}^{N_n} \hat{\rho}_a(x)$, where
$\hat{\rho}_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 in the system.
Under the BO approximation, $\hat{\rho}_a(x) = Z_a e\, \delta^{(3)}(\vec{r} - \vec{R}_a)$,
where $Z_a$ is the nucleus $a$'s atomic number, and $\vec{R}_a$
should be understood as the direct sum of the position of each nucleus of type $a$ (see Appendix of Ref.~\cite{Ichikawa2013} 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 $\hat{\vec{j}}(x)=\hat{\vec{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 \ref{sec:H}.

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
\cite{Jackson,Loudon,LLCTF},
and express $\hat{A}_\mu(x)$ by such solutions \cite{Tachibana2003,Ichikawa2013}.
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 \cite{Yang1950}.
We adopt the Coulomb gauge, $\vec{\nabla} \cdot \hat{\vec{A}}(x) = 0$.
Then, the scalar potential is given by
\begin{eqnarray}
\hat{A}_0(ct, \vec{r}) = \int d^3\vec{s}\, \frac{\hat{\rho}(ct, \vec{s})}{|\vec{r}-\vec{s}|}, \label{eq:A0}
\end{eqnarray}
and the vector potential by
\begin{eqnarray}
\hat{\vec{A}}(ct, \vec{r}) = \hat{\vec{A}}_{\rm rad}(ct, \vec{r}) + \hat{\vec{A}}_A(ct, \vec{r}), \label{eq:vecA}
\end{eqnarray}
where the first term is the quantized free radiation field in the Coulomb gauge \cite{Loudon,SakuraiQFT} and the second term is the retarded potential.
$\hat{\vec{A}}_{\rm rad}(ct, \vec{r})$ is expressed as
\begin{eqnarray}
\hat{A}^k_{\rm rad}(ct, \vec{r})
 &=& \frac{\sqrt{4\pi \hbar^2 c}}{\sqrt{(2\pi\hbar)^3}} \sum_{\sigma=\pm1} \int \frac{d^3 \vec{p}}{\sqrt{2p^0}}
\bigg[ \hat{a}_{\vec{p}_\sigma} e^k(\vec{p},\sigma) e^{-i c p^0 t/\hbar} e^{i \vec{p}\cdot \vec{r} /\hbar} \nonumber \\
& &+ \hat{a}^\dagger_{\vec{p}_\sigma} e^{* k}(\vec{p},\sigma) e^{i c p^0 t/\hbar} e^{-i \vec{p}\cdot \vec{r} /\hbar} \bigg] ,
\label{eq:Aradvec}
\end{eqnarray}
where $\hat{a}_{\vec{p}_\sigma}$ ($\hat{a}^\dagger_{\vec{p}_\sigma}$) are the annihilation (creation) operator
of the photon with momentum $\vec{p}$ and helicity $\sigma$, and
$\vec{e}$ is the polarization vector. They satisfy the commutation relation
$[ \hat{a}_{\vec{p}_\sigma}, \hat{a}^\dagger_{\vec{q}_\tau} ] = \delta(\vec{p}-\vec{q}) \delta_{\sigma \tau}$,
and commutators of the other combinations are zero.
Our convention for $\vec{e}(\vec{p},\sigma)$ is found in Ref.~\cite{Ichikawa2013}.
$\hat{\vec{A}}_A(ct, \vec{r})$ is expressed as
\begin{eqnarray}
\hat{\vec{A}}_A(ct, \vec{r}) = \frac{1}{c} \int d^3\vec{s}\, \frac{\hat{\vec{j}}_T(c u, \vec{s})}{|\vec{r}-\vec{s}|},  \label{eq:A_A}
\end{eqnarray}
where we define the retarded time
$u = t - |\vec{r}-\vec{s}|/c$, and $\hat{\vec{j}}_T(x)$ is the transversal component
of the current.
Writing explicitly,
\begin{eqnarray}
\hat{j}_T^k(cu,\vec{s}) &=&  \sum_{p,q=1}^{N_D} \sum_{c,d=\pm} \left\{
 j^k_{p^c q^d}(\vec{s}) \hat{\cal E}_{p^c q^d}(u)
+ E^k_{p^c q^d}(\vec{s}) \frac{d \hat{\cal E}_{p^c q^d}}{dt}(u)
\right\},   \label{eq:jT}
\end{eqnarray}
where
\begin{eqnarray}
E^{k}_{n^a m^b}(\vec{R})
=
-\frac{Z_e e}{4\pi} \int d^3\vec{s}\, \psi_{n^a}^\dagger(\vec{s}) \psi_{m^b}(\vec{s})  \frac{(\vec{s}-\vec{R})^k}{|\vec{s}-\vec{R}|^3}.   \label{eq:ef_int}
\end{eqnarray}

As we mentioned in Sec.~\ref{sec:introduction}, we do not consider the usual asymptotic states of the field operator,
in which the state becomes non-interacting as $t \rightarrow -\infty$ asymptotically.
We assume that the interaction is absent for $t<0$ and the free part expressed above by $\hat{A}_{\rm rad, \mu}(x)$
is considered to realize at $t = 0$.
This, in particular, implies the initial condition
\begin{eqnarray}
\hat{j}^\mu(cu,\vec{s}) = 0, \quad u<0,
\end{eqnarray}
and enables us to rewrite the retarded potential as \cite{Tachibana2013}
\begin{eqnarray}
\hat{\vec{A}}_A(ct, \vec{r}) &=&
\frac{1}{c^2 \pi} \int_{0}^{t }du  \int_{-\infty}^{\infty}d\alpha \exp \left( i\alpha (t-u)^2 \right)
\int d^3\vec{s}\, \hat{\vec{j}}_T(c u, \vec{s}) \exp \left( -i\alpha  \frac{(\vec{r}-\vec{s})^2}{c^2} \right).
\label{eq:A_A_mod}
\end{eqnarray}
To derive this, we use the following formulae for the delta function
$\delta(x^2-a^2) = \left\{ \delta(x-a) + \delta(x+a) \right\}/(2a)$ with $a>0$,
and
$\delta( (t-u)^2 - (\vec{r}-\vec{s})^2/c^2 ) = \frac{1}{2\pi} \int_{-\infty}^{\infty} d\alpha \exp [ i\alpha \{ (t-u)^2 -(\vec{r}-\vec{s})^2/c^2\} ]$.
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.~\cite{Tachibana2005}, 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
$m_e = m_{eB} + \delta m_e$, $e = \sqrt{Z_3} e_B$,
$\hat{\psi}(x) = \hat{\psi}_B(x)/ \sqrt{Z_2}$ and $\hat{A}^\mu(x) = \hat{A}_B^\mu(x)/ \sqrt{Z_3}$,
where ``B" in the subscripts denotes a bare parameter or a bare field,
and $ \delta 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.

\subsection{Time evolution of annihilation and creation operators}  \label{sec:qnum_evac}
The time derivative of $\hat{e}_{n^a}$ is given by substituting the expansion Eq.~\eqref{eq:psi_expand2}
into the Dirac field equation, multiplying by $\psi^\dagger_{n^a}(\vec{r})$, integrating over $\vec{r}$,
and using the orthonormality condition \cite{Ichikawa2013}.
This leads to
\begin{eqnarray}
\hspace{-0.6cm}
i\hbar \frac{\partial \hat{e}_{n^a}}{\partial t} &=&
  \sum_{m=1}^{N_D}  \sum_{b=\pm}  h_{n^a m^b} \hat{e}_{m^b}
  +    \sum_{m, p, q=1}^{N_D} \sum_{b, c, d=\pm}(n^a m^b | p^c q^d) \hat{e}^\dagger_{p^c} \hat{e}_{q^d}  \hat{e}_{m^b}
  \nonumber \\
&-&\frac{1}{c^2} \sum_{m=1}^{N_D}  \sum_{b=\pm}  \sum_{k=1}^3 \int d^3\vec{r}\, d^3\vec{s}\,
j^k_{n^a m^b}(\vec{r})  \frac{\hat{j}^k_T(c u, \vec{s})}{|\vec{r}-\vec{s}|}  \hat{e}_{m^b}
 \nonumber \\
&-& \frac{\sqrt{4\pi \hbar^2}}{\sqrt{c (2\pi\hbar)^3}}
 \sum_{m=1}^{N_D}  \sum_{b=\pm} \sum_{k=1}^3  \sum_{\sigma=\pm1} \int \frac{d^3 \vec{p}}{\sqrt{2p^0}}
\bigg[
F^k _{n^a m^b}(\vec{p}) e^k(\vec{p},\sigma) e^{-i c p^0 t/\hbar}\, \hat{a}_{\vec{p}_\sigma} \hat{e}_{m^b} \nonumber \\
& &+F^k_{n^a m^b}(-\vec{p}) e^{* k}(\vec{p},\sigma) e^{i c p^0 t/\hbar}\,   \hat{a}^\dagger_{\vec{p}_\sigma}   \hat{e}_{m^b}
\bigg],
\label{eq:dedt_BO}
\end{eqnarray}
where the various coefficient matrices are defined as follows.
In the first term, we define
\begin{eqnarray}
h_{n^a m^b} = T_{n^a m^b} + M_{n^a m^b} +   \sum_{a=1}^{N_n} (Z_a e) V_{n^a m^b}(\vec{R}_a),
\label{eq:inth}
\end{eqnarray}
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}(\vec{R})$ is the nuclear attraction integral, respectively defined as
\begin{eqnarray}
T_{n^a m^b}
&\equiv&  -i\hbar c \int d^3\vec{r}\,\psi_{n^a}^\dagger(\vec{r}) \gamma^0 \vec{\gamma} \cdot \vec{\nabla}  \psi_{m^b}(\vec{r}), \label{eq:kin_int} \\
M_{n^a m^b}
&\equiv& m_e c^2 \int d^3\vec{r}\,\psi_{n^a}^\dagger(\vec{r}) \gamma^0  \psi_{m^b}(\vec{r}), \label{eq:mass_int} \\
V_{n^a m^b}(\vec{R})
&\equiv& (Z_e  e) \int d^3\vec{s}\, \frac{\psi_{n^a}^\dagger(\vec{s}) \psi_{m^b}(\vec{s}) }{|\vec{s}-\vec{R}|}.   \label{eq:nuc_int}
\end{eqnarray}
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.~\eqref{eq:dedt_BO}, we define the electronic repulsion integral
\begin{eqnarray}
(n^a m^b | p^c q^d)
\equiv
(Z_e e)^2 \int d^3\vec{r}\, d^3\vec{s}\, \psi_{n^a}^\dagger(\vec{r}) \psi_{m^b}(\vec{r})
\frac{1}{|\vec{r}-\vec{s}|} \psi_{p^c}^\dagger(\vec{s}) \psi_{q^d}(\vec{s}).
   \label{eq:twoele_int}
\end{eqnarray}
This term originates in the electronic charge contribution to the scalar potential term.
The third and forth terms of Eq.~\eqref{eq:dedt_BO} originate in the vector potential term.
The former comes from the retarded part $\hat{\vec{A}}_A$ \eqref{eq:A_A}
and the latter from the free radiation part $\hat{\vec{A}}_{\rm rad}$ \eqref{eq:Aradvec}.
Finally, in the forth term of Eq.~\eqref{eq:dedt_BO}, we define
\begin{eqnarray}
F^k_{n^a m^b}(\vec{p}) \equiv \int d^3\vec{r}\,  j^k_{n^a m^b}(\vec{r}) e^{i \vec{p}\cdot \vec{r} /\hbar}.
\end{eqnarray}

\subsection{Time evolution of excitation and $\hat{e}^\dagger \hat{a} \hat{e}$-type operators}  \label{sec:qnum_excitation}
In this paper, we are interested in the time evolution of electronic charge density, and its operator $\hat{\rho}_e(x)$ is
expressed by the excitation operator as Eq.~\eqref{eq:op_rhoe}.
Since the excitation operator carries every time dependence of $\hat{\rho}_e(x)$,
what we need to know is the time evolution equation of the excitation operator.

The time derivative of the excitation operator $\hat{\cal E}_{n^a m^b}$ \eqref{eq:op_excitation} can be written as
\begin{eqnarray}
\frac{\partial  \hat{\cal E}_{n^a m^b}}{\partial t}
=
 (\hat{O}_{m^b n^a})^\dagger + \hat{O}_{n^a m^b},  \label{eq:dcalEdt}
\end{eqnarray}
where we define
\begin{eqnarray}
\hat{O}_{n^a m^b}  \equiv \hat{e}^\dagger_{n^a}   \frac{\partial \hat{e}_{m^b}}{\partial t}. \label{eq:def_op_O}
\end{eqnarray}
Using Eq.~\eqref{eq:dedt_BO}, this can be readily obtained as
\begin{eqnarray}
i\hbar \hat{O}_{n^a m^b}
&=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} \hat{e}^\dagger_{n^a} \hat{e}_{r^e}
  +    \sum_{r, p, q=1}^{N_D} \sum_{e, c, d=\pm}(m^b r^e | p^c q^d) \hat{e}^\dagger_{n^a}   \hat{e}^\dagger_{p^c} \hat{e}_{q^d}    \hat{e}_{r^e}
  \nonumber \\
&-&\frac{1}{c^2} \sum_{r=1}^{N_D}  \sum_{e=\pm}  \sum_{k=1}^3 \int d^3\vec{r}\, d^3\vec{s}\,
j^k_{m^b r^e}(\vec{r}) \hat{e}^\dagger_{n^a} \frac{\hat{j}^k_T(c u, \vec{s})}{|\vec{r}-\vec{s}|}  \hat{e}_{r^e}
 \nonumber \\
&-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm} \sum_{\sigma=\pm1} \int \frac{d^3 \vec{p}}{\sqrt{2p^0}}
 \left[ {\cal F}_{m^b r^e \vec{p}_\sigma}(t) \hat{e}^\dagger_{n^a}\hat{a}_{\vec{p}_\sigma} \hat{e}_{r^e}
+{\cal F}^*_{r^e m^b \vec{p}_\sigma}(t) \hat{e}^\dagger_{n^a} \hat{a}^\dagger_{\vec{p}_\sigma}   \hat{e}_{r^e}
\right],
\nonumber \\
\label{eq:op_O}
\end{eqnarray}
where we have defined
\begin{eqnarray}
{\cal F}_{n^a m^b \vec{p}_\sigma}(t) \equiv  \sum_{k=1}^3 F^k _{n^a m^b}(\vec{p}) e^k(\vec{p},\sigma) e^{-i c p^0 t/\hbar},
\end{eqnarray}
to make the expression shorter.
For later convenience, we call the second term of Eq.~\eqref{eq:op_O} ``four-electron term", the third term ``retarded potential term",
and the forth term ``radiation term".

Next, we consider the time derivative of the $\hat{e}^\dagger \hat{a} \hat{e}$-type operators.
Since the time derivative of $\hat{e}^\dagger_{m^b}\hat{a}^\dagger_{\vec{p}_\sigma} \hat{e}_{n^a}$ is known by
taking the Hermite conjugate of that of $\hat{e}^\dagger_{n^a}\hat{a}_{\vec{p}_\sigma} \hat{e}_{m^b}$,
we only show the latter below.
This can be expressed as
\begin{eqnarray}
\frac{\partial }{\partial t} \left\{ \hat{e}^\dagger_{n^a}\hat{a}_{\vec{p}_\sigma}\hat{e}_{m^b} \right\}
=
(\hat{Q}_{m^b \vec{p}_\sigma n^a})^\dagger + \hat{P}_{n^a \vec{p}_\sigma m^b}, \label{eq:op_deaedt}
\end{eqnarray}
where we define
\begin{eqnarray}
 \hat{P}_{n^a \vec{p}_\sigma m^b} &\equiv& \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma} \frac{\partial\hat{e}_{m^b}}{\partial t},  \label{eq:def_opP} \\
  \hat{Q}_{n^a \vec{p}_\sigma m^b} &\equiv& \hat{e}^\dagger_{n^a} \hat{a}^\dagger_{\vec{p}_\sigma} \frac{\partial\hat{e}_{m^b}}{\partial t}. \label{eq:def_opQ}
\end{eqnarray}
Using Eq.~\eqref{eq:dedt_BO}, we obtain
\begin{eqnarray}
i\hbar \hat{P}_{n^a \vec{p}_\sigma m^b}
&=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma} \hat{e}_{r^e}
  +    \sum_{r, p, q=1}^{N_D} \sum_{e, c, d=\pm}(m^b r^e | p^c q^d) \hat{e}^\dagger_{n^a}\hat{a}_{\vec{p}_\sigma}  \hat{e}^\dagger_{p^c} \hat{e}_{q^d}    \hat{e}_{r^e}
  \nonumber \\
&-&\frac{1}{c^2} \sum_{r=1}^{N_D}  \sum_{e=\pm}  \sum_{k=1}^3 \int d^3\vec{r}\, d^3\vec{s}\,
j^k_{m^b r^e}(\vec{r}) \hat{e}^\dagger_{n^a}\hat{a}_{\vec{p}_\sigma} \frac{\hat{j}^k_T(c u, \vec{s})}{|\vec{r}-\vec{s}|}  \hat{e}_{r^e}
 \nonumber \\
&-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm} \sum_{\tau=\pm1} \int \frac{d^3 \vec{q}}{\sqrt{2q^0}} \nonumber \\
&& \times  \left[ {\cal F}_{m^b r^e \vec{q}_\tau}(t)  \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{a}_{\vec{q}_\tau} \hat{e}_{r^e}
+{\cal F}^*_{r^e m^b \vec{q}_\tau}(t)   \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{a}^\dagger_{\vec{q}_\tau}   \hat{e}_{r^e}
\right],  \label{eq:op_P}
\end{eqnarray}
and
\begin{eqnarray}
i\hbar \hat{Q}_{n^a \vec{p}_\sigma m^b}
&=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} \hat{e}^\dagger_{n^a} \hat{a}^\dagger_{\vec{p}_\sigma}  \hat{e}_{r^e}
  +    \sum_{r, p, q=1}^{N_D} \sum_{e, c, d=\pm}(m^b r^e | p^c q^d) \hat{e}^\dagger_{n^a}\hat{a}^\dagger_{\vec{p}_\sigma}    \hat{e}^\dagger_{p^c} \hat{e}_{q^d}    \hat{e}_{r^e}
  \nonumber \\
&-&\frac{1}{c^2} \sum_{r=1}^{N_D}  \sum_{e=\pm}  \sum_{k=1}^3 \int d^3\vec{r}\, d^3\vec{s}\,
j^k_{m^b r^e}(\vec{r}) \hat{e}^\dagger_{n^a}\hat{a}^\dagger_{\vec{p}_\sigma}   \frac{\hat{j}^k_T(c u, \vec{s})}{|\vec{r}-\vec{s}|}  \hat{e}_{r^e}
 \nonumber \\
&-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm} \sum_{\tau=\pm1} \int \frac{d^3 \vec{q}}{\sqrt{2q^0}} \nonumber \\
& & \times \left[ {\cal F}_{m^b r^e \vec{q}_\tau}(t)  \hat{e}^\dagger_{n^a} \hat{a}^\dagger_{\vec{p}_\sigma}   \hat{a}_{\vec{q}_\tau}  \hat{e}_{r^e}
 +{\cal F}^*_{r^e m^b \vec{q}_\tau}(t) \hat{e}^\dagger_{n^a} \hat{a}^\dagger_{\vec{p}_\sigma}   \hat{a}^\dagger_{\vec{q}_\tau}    \hat{e}_{r^e}
\right].\label{eq:op_Q}
\end{eqnarray}
We note that the last term of $\hat{P}_{n^a \vec{p}_\sigma m^b}$ includes the operator product of the form
$\hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{a}^\dagger_{\vec{q}_\tau}   \hat{e}_{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.~\ref{sec:selfenergy}.

\section{Evolution equations for density matrix}  \label{sec:dm}

\subsection{Definition of density matrix} \label{sec:def_dm}
We begin by introducing notations regarding expectation values.
We denote the expectation value of the excitation operator with respect to the Heisenberg initial ket $| \Phi \rangle$ by
\begin{eqnarray}
{\cal E}_{n^a m^b} \equiv  \langle \Phi |   \hat{\cal E}_{n^a m^b}   | \Phi \rangle,
\end{eqnarray}
and call this quantity the density matrix.
Below, we sometimes write just brackets around the operator to denote the expectation value with respect to $| \Phi \rangle$,
namely, $\langle  \cdots \rangle \equiv  \langle \Phi | \cdots | \Phi \rangle $ where $\cdots$ stands for some operators.
Since $(\hat{\cal E}_{n^a m^b})^\dagger = \hat{e}^\dagger_{m^b} \hat{e}_{n^a}  = \hat{\cal E}_{m^b n^a}$,
taking the expectation value yields $({\cal E}_{n^a m^b})^* = {\cal E}_{m^b n^a}$,
showing that the density matrix is a Hermite matrix.
Similarly, we define the expectation value of Eq.~\eqref{eq:def_op_O} as
${\cal O}_{n^a m^b} \equiv   \langle \Phi |  \hat{O}_{n^a m^b}  | \Phi \rangle$.
Then, since
 $\langle  (\hat{O}_{n^a m^b})^\dagger \rangle = ({\cal O}_{n^a m^b})^* = {\cal O}^\dagger_{m^b n^a}$,
taking the expectation value of Eq.~\eqref{eq:dcalEdt} leads to
\begin{eqnarray}
\frac{\partial  {\cal E}_{n^a m^b}}{\partial t} =  {\cal O}^\dagger_{n^a m^b} + {\cal O}_{n^a m^b},   \label{eq:ev_dcalEdt}
\end{eqnarray}
which is the time evolution equation for the density matrix.

Using the density matrix, the expectation value of the electronic charge density operator, Eq.~\eqref{eq:op_rhoe},
can be written as
\begin{eqnarray}
 \langle \Phi |  \hat{\rho}_e(x)| \Phi \rangle =  \sum_{n,m=1}^{N_D} \sum_{a,b=\pm} \rho_{n^a m^b}(\vec{r}) {\cal E}_{n^a m^b}(t).
 \end{eqnarray}
However, this quantity has a non-zero value when the Heisenberg initial ket
is the vacuum, $| \Phi \rangle = | 0 \rangle$, at $t=0$.
(Remember that $\hat{e}_{n^-} $ is the positron creation operator.)
As is done in the ordinary QED, this can be remedied by computing the
expectation value of the normal-ordered product of the operators at $t=0$.
In the case of $\hat{\rho}_e(x)$, this can be accomplished by subtracting the vacuum expectation value
of $\hat{\rho}_e(x)$ at $t=0$.
Then, we shall define the electronic charge density $\rho_e(x)$ as the expectation value
of the electronic charge density operator after this subtraction.
Writing explicitly,
\begin{eqnarray}
\rho_e(x) = \sum_{n,m=1}^{N_D} \sum_{a,b=\pm} \rho_{n^a m^b}(\vec{r})
\left\{ {\cal E}_{n^a m^b}(t) - {\cal E}^0_{n^a m^b}(t=0) \right\},   \label{eq:phys_rhoe}
\end{eqnarray}
where we define ${\cal E}^0_{n^a m^b} \equiv  \langle 0|   \hat{\cal E}_{n^a m^b}   | 0 \rangle$.

\subsection{Four-electron term and retarded potential term} \label{sec:approx_4e_ret}
The four-electron term consists of four electron creation or annihilation operators.
We approximate this term by
\begin{eqnarray}
\langle \hat{e}^\dagger_{n^a}  \hat{e}^\dagger_{p^c}    \hat{e}_{q^d}  \hat{e}_{r^e}   \rangle
&\approx&
\langle \hat{e}^\dagger_{n^a}  \hat{e}_{r^e} \rangle  \langle \hat{e}^\dagger_{p^c}    \hat{e}_{q^d}    \rangle
- \langle \hat{e}^\dagger_{n^a}  \hat{e}_{q^d}  \rangle  \langle \hat{e}^\dagger_{p^c}     \hat{e}_{r^e}   \rangle \nonumber \\
&=&
{\cal E}_{n^a r^e} {\cal E}_{p^c q^d} -{\cal E}_{n^a q^d} {\cal E}_{p^c r^e}.
\label{eq:approx_4e}
\end{eqnarray}
This decomposition holds exactly for $t=0$ and the approximation is motivated by this fact.
In our previous paper Ref.~\cite{Ichikawa2013}, 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 $\hat{j}^k_T$ is expressed as Eq.~\eqref{eq:jT}, it seems to
have the same structure as the four-electron term.
However, since $\hat{j}^k_T$ is computed at the retarded time and contains a time derivative of
the excitation operator, the approximation like Eq.~\eqref{eq:approx_4e} is not applicable.
Thus, we approximate this term by just replacing $\hat{j}^k_T$ by its expectation value as
\begin{eqnarray}
\langle \hat{e}^\dagger_{n^a}(t) \hat{j}^k_T(c u, \vec{s}) \hat{e}_{r^e}(t) \rangle
\approx
\langle  \hat{j}^k_T(c u, \vec{s})\rangle \langle   \hat{e}^\dagger_{n^a}(t)\hat{e}_{r^e}(t) \rangle
=
{\cal J}^k_T(c u, \vec{s}) {\cal E}_{n^a r^e}(t),
\label{eq:approx_jT}
\end{eqnarray}
where we define ${\cal J}^k_T \equiv   \langle \Phi | \hat{j}^k_T  | \Phi \rangle$.
With this approximation,
the two space integrations over $\vec{r}$ and $\vec{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 $\vec{r}$ and $\vec{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.~\eqref{eq:A_A_mod}.
Then, for the expectation value of the retarded term, we have
\begin{eqnarray}
-\frac{1}{c^2} \sum_{r=1}^{N_D}  \sum_{e=\pm}  \sum_{k=1}^3 \int d^3\vec{r}\, d^3\vec{s}\,
j^k_{m^b r^e}(\vec{r}) \frac{{\cal J}^k_T(c u, \vec{s})}{|\vec{r}-\vec{s}|}   {\cal E}_{n^a r^e}(t)
= \sum_{r=1}^{N_D}  \sum_{e=\pm}  I_{j_T}[{\cal E},\dot{\cal E}]_{m^b r^e}(t)  {\cal E}_{n^a r^e}(t), \label{eq:approx_O3rd}
\end{eqnarray}
where we define
\begin{eqnarray}
I_{j_T}[{\cal E},\dot{\cal E}]_{m^b r^e}(t) &\equiv&
-\frac{1}{c^3\pi} \sum_{p,q=1}^{N_D}  \sum_{c,d=\pm} \int_{0}^{t }du  \int_{-\infty}^{\infty}d\alpha \exp \left( i\alpha (t-u)^2 \right)  \nonumber \\
& &\times \left\{
I_{jj, m^b r^e p^c q^d}(\alpha) {\cal E}_{p^c q^d}(u)  + I_{jE, m^b r^e p^c q^d}(\alpha)  \frac{d {\cal E}_{p^c q^d}}{dt}(u)
\right\},  \label{eq:def_IjT}
\end{eqnarray}
with
\begin{eqnarray}
I_{jj, m^b r^e p^c q^d}(\alpha) &\equiv&
\sum_{k=1}^3  \int d^3\vec{r}\,  d^3\vec{s}\,  j^k_{m^b r^e}(\vec{r})  j^k_{p^c q^d}(\vec{s}) \exp \left( -i\alpha  \frac{(\vec{r}-\vec{s})^2}{c^2} \right),
\label{eq:def_Ijj} \\
I_{jE, m^b r^e p^c q^d}(\alpha) &\equiv&
\sum_{k=1}^3  \int d^3\vec{r}\,  d^3\vec{s}\,  j^k_{m^b r^e}(\vec{r}) E^k_{p^c q^d}(\vec{s}) \exp \left( -i\alpha  \frac{(\vec{r}-\vec{s})^2}{c^2} \right).
\label{eq:def_IjE}
\end{eqnarray}
Here, we put $[{\cal E},\dot{\cal E}]$ after $I_{j_T}$ for the notation of the integral \eqref{eq:def_IjT}
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.~\eqref{eq:def_Ijj} and \eqref{eq:def_IjE} are defined by Eqs.~\eqref{eq:jk} and \eqref{eq:ef_int},
$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~\ref{sec:molint}.

\subsection{Radiation term} \label{sec:approx_rad}
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.~\cite{Ichikawa2013}.
In this approximation, we use
\begin{eqnarray}
\langle \hat{e}^\dagger_{n^a}\hat{a}_{\vec{p}_\sigma}  \hat{e}_{r^e} \rangle
&\approx&
 \langle \hat{a}_{\vec{p}_\sigma}  \rangle  \langle \hat{e}^\dagger_{n^a} \hat{e}_{r^e} \rangle
= \langle \hat{a}_{\vec{p}_\sigma}  \rangle  {\cal E}_{n^a r^e},
\end{eqnarray}
and, similarly, $\langle \hat{e}^\dagger_{n^a}\hat{a}^\dagger_{\vec{p}_\sigma} \hat{e}_{r^e} \rangle
\approx \langle \hat{a}_{\vec{p}_\sigma}  \rangle ^* {\cal E}_{n^a r^e}$.
This factorization holds exactly for $t=0$, when $ \hat{e}_{n^a}$ and $\hat{a}_{\vec{p}_\sigma} $ 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
\begin{eqnarray}
i\hbar {\cal O}_{n^a m^b}
&=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} {\cal E}_{n^a r^e}
  +    \sum_{r, p, q=1}^{N_D} \sum_{e, c, d=\pm}(m^b r^e | p^c q^d) \left( {\cal E}_{n^a r^e} {\cal E}_{p^c q^d} -{\cal E}_{n^a q^d} {\cal E}_{p^c r^e} \right)
  \nonumber \\
&+& \sum_{r=1}^{N_D}  \sum_{e=\pm}  I_{j_T}[{\cal E},\dot{\cal E}]_{m^b r^e}(t)  {\cal E}_{n^a r^e}
 \nonumber \\
&-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm} \sum_{\sigma=\pm1} \int \frac{d^3 \vec{p}}{\sqrt{2p^0}}
\left[ {\cal F}_{m^b r^e \vec{p}_\sigma}(t) \langle \hat{a}_{\vec{p}_\sigma} \rangle  {\cal E}_{n^a r^e}
+{\cal F}^*_{r^e m^b \vec{p}_\sigma}(t)\langle \hat{a}_{\vec{p}_\sigma} \rangle ^* {\cal E}_{n^a r^e}
\right]. \label{eq:calO_1} \nonumber \\
\end{eqnarray}
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 $\hat{e}^\dagger \hat{a} \hat{e}$-type operators
simultaneously with one for the density matrix.
We first define
\begin{eqnarray}
{\cal E}_{n^a \vec{p}_\sigma m^b} \equiv  \langle \Phi |   \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{e}_{m^b}    | \Phi \rangle.
\end{eqnarray}
Note that $\langle \Phi |   \hat{e}^\dagger_{n^a} \hat{a}^\dagger_{\vec{p}_\sigma} \hat{e}_{m^b} | \Phi \rangle = {\cal E}^*_{m^b \vec{p}_\sigma n^a}$.
We next define the expectation value of the operator $\hat{P}_{n^a \vec{p}_\sigma m^b}$, Eq.~\eqref{eq:def_opP},
as ${\cal P}_{n^a \vec{p}_\sigma m^b} \equiv   \langle \Phi |  \hat{P}_{n^a \vec{p}_\sigma m^b}  | \Phi \rangle$.
Similarly for the operator $\hat{Q}_{n^a \vec{p}_\sigma m^b}$, Eq.~\eqref{eq:def_opQ}, we define ${\cal Q}_{n^a \vec{p}_\sigma m^b} \equiv   \langle \Phi |  \hat{Q}_{n^a \vec{p}_\sigma m^b}  | \Phi \rangle$.
Then, the expectation value of Eq.~\eqref{eq:op_O} can be written as
\begin{eqnarray}
i\hbar {\cal O}_{n^a m^b}
&=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} {\cal E}_{n^a r^e}
  +    \sum_{r, p, q=1}^{N_D} \sum_{e, c, d=\pm}(m^b r^e | p^c q^d) \left( {\cal E}_{n^a r^e} {\cal E}_{p^c q^d} -{\cal E}_{n^a q^d} {\cal E}_{p^c r^e} \right)
  \nonumber \\
&+& \sum_{r=1}^{N_D}  \sum_{e=\pm}  I_{j_T}[{\cal E},\dot{\cal E}]_{m^b r^e}(t)  {\cal E}_{n^a r^e}
 \nonumber \\
&-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm}  \sum_{\sigma=\pm1} \int \frac{d^3 \vec{p}}{\sqrt{2p^0}}
\left[ {\cal F}_{m^b r^e \vec{p}_\sigma}(t) {\cal E}_{n^a \vec{p}_\sigma r^e}
+{\cal F}^*_{r^e m^b \vec{p}_\sigma}(t)  {\cal E}^*_{r^e \vec{p}_\sigma n^a}
\right], \label{eq:calO_2} \nonumber \\
\end{eqnarray}
and the time evolution of  ${\cal E}_{n^a \vec{p}_\sigma m^b}$ can be expressed as
\begin{eqnarray}
\frac{\partial {\cal E}_{n^a \vec{p}_\sigma m^b}}{\partial t}
=
{\cal Q}^*_{m^b \vec{p}_\sigma n^a} + {\cal P}_{n^a \vec{p}_\sigma m^b}, \label{eq:ev_deaedt}
\end{eqnarray}
which is obtained by taking the expectation value of Eq.~\eqref{eq:op_deaedt}.

As for ${\cal P}_{n^a \vec{p}_\sigma m^b}$, the expectation value of the operator product in the second term of Eq.~\eqref{eq:op_P} is approximated to be
\begin{eqnarray}
\langle \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{e}^\dagger_{p^c} \hat{e}_{q^d}    \hat{e}_{r^e}  \rangle
\approx
\langle  \hat{a}_{\vec{p}_\sigma} \rangle \left\{
\langle \hat{e}^\dagger_{n^a}  \hat{e}_{r^e} \rangle  \langle \hat{e}^\dagger_{p^c}    \hat{e}_{q^d}    \rangle
- \langle \hat{e}^\dagger_{n^a}  \hat{e}_{q^d}  \rangle  \langle \hat{e}^\dagger_{p^c}     \hat{e}_{r^e}   \rangle
\right\}.
\end{eqnarray}
In this approximation, we assume that the initial photon state is a number state and not a coherent state.
Then, since $\langle  \hat{a}_{\vec{p}_\sigma} \rangle=0$, the contribution of this term vanishes.
The third term of Eq.~\eqref{eq:op_P} includes the retarded time and its expectation value would be
approximated in the same way as the retarded potential term of Eq.~\eqref{eq:op_O}.
Following the procedure described below Eq.~\eqref{eq:approx_jT},
we have the expression of the form of Eq.~\eqref{eq:approx_O3rd} with ${\cal E}_{n^a r^e}(t)$
replaced by ${\cal E}_{n^a \vec{p}_\sigma r^e}(t)$.
The forth term of Eq.~\eqref{eq:op_P} includes two four-operator terms, which are approximated to be
\begin{eqnarray}
\langle \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{a}_{\vec{q}_\tau} \hat{e}_{r^e} \rangle
 &\approx&
\langle  \hat{a}_{\vec{p}_\sigma}  \hat{a}_{\vec{q}_\tau}  \rangle \langle \hat{e}^\dagger_{n^a} \hat{e}_{r^e} \rangle,
 \\
\langle \hat{e}^\dagger_{n^a} \hat{a}_{\vec{p}_\sigma}  \hat{a}^\dagger_{\vec{q}_\tau}   \hat{e}_{r^e} \rangle
 &\approx&
 \langle  \hat{a}_{\vec{p}_\sigma}  \hat{a}^\dagger_{\vec{q}_\tau}  \rangle  \langle \hat{e}^\dagger_{n^a} \hat{e}_{r^e} \rangle.
\end{eqnarray}
Since we do not consider a coherent state for the initial photon state as mentioned above,
$\langle  \hat{a}_{\vec{p}_\sigma}  \hat{a}_{\vec{q}_\tau}  \rangle=0$  in the first equation
and its contribution vanishes.
In the second equation, since
$\langle  \hat{a}_{\vec{p}_\sigma}  \hat{a}^\dagger_{\vec{q}_\tau}  \rangle
= (n_{\vec{p}_\sigma}+1) \delta^{(3)}(\vec{p}-\vec{q}) \delta_{\sigma \tau}$,
where $n_{\vec{p}_\sigma}$ is the occupation number of the photon mode $(\vec{p},\sigma)$
in the initial state, it may give a non-zero contribution.
Putting these approximations together,
we have
\begin{eqnarray}
i\hbar {\cal P}_{n^a \vec{p}_\sigma m^b} &=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} {\cal E}_{n^a  \vec{p}_\sigma r^e}
  + \sum_{r=1}^{N_D}  \sum_{e=\pm}  I_{j_T}[{\cal E},\dot{\cal E}]_{m^b r^e}(t)  {\cal E}_{n^a \vec{p}_\sigma  r^e}
 \nonumber \\
 &-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm} \frac{1}{\sqrt{2p^0}} {\cal F}^*_{r^e m^b \vec{p}_\sigma}(t) (n_{\vec{p}_\sigma}+1) {\cal E}_{n^a  r^e},
 \label{eq:calP}
\end{eqnarray}
when we assume a number state for the photon initial state.
The expectation value of Eq.~\eqref{eq:op_Q} can be approximated in a similar manner to be
\begin{eqnarray}
i\hbar {\cal Q}_{n^a \vec{p}_\sigma m^b} &=&
  \sum_{r=1}^{N_D}  \sum_{e=\pm}  h_{m^b r^e} {\cal E}^*_{r^e  \vec{p}_\sigma n^a}
  + \sum_{r=1}^{N_D}  \sum_{e=\pm}  I_{j_T}[{\cal E},\dot{\cal E}]_{m^b r^e}(t) {\cal E}^*_{r^e  \vec{p}_\sigma n^a}
 \nonumber \\
 &-& \frac{1}{\sqrt{2\pi^2 \hbar c}}
 \sum_{r=1}^{N_D}  \sum_{e=\pm} \frac{1}{\sqrt{2p^0}} {\cal F}_{m^b r^e \vec{p}_\sigma}(t) n_{\vec{p}_\sigma} {\cal E}_{n^a  r^e},
 \label{eq:calQ}
\end{eqnarray}
where we have used
$\langle  \hat{a}^\dagger_{\vec{p}_\sigma}  \hat{a}_{\vec{q}_\tau}  \rangle
= n_{\vec{p}_\sigma} \delta^{(3)}(\vec{p}-\vec{q}) \delta_{\sigma \tau}$.

\section{Results}  \label{sec: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 \cite{QEDynamics} developed in our group.
In Sec.~\ref{sec:calc_setup}, we describe our setups for numerical calculation
including the initial condition for the density matrix.
The results of two approximation methods discussed in Sec.~\ref{sec:approx_rad}
are respectively presented in Sec.~\ref{sec:coherent} and Sec.~\ref{sec:selfenergy}.

\subsection{Setups for numerical calculation} \label{sec:calc_setup}
To perform  numerical calculation, we first need to determine an orthonormal set of
expansion functions to define the electron creation and annihilation operators,
as explained in Sec.~\ref{sec:qnum_defac}.
We generate the set by solving the Dirac equation with the four-component Dirac-Coulomb hamiltonian.
They are computed by the publicly available program package DIRAC \cite{DIRAC12},
using the Hartree-Fock method with the STO-3G basis set.
Note that, for the hydrogen atom in this basis set, we have two ($=N_D$) orbitals for
electron and positron respectively taking into account the Kramers partners.
The density matrix is 4 $\times$ 4 matrix whose components denote
electron ($1^+$), its Kramers partner ($\bar{1}^+$), positron ($1^-$), and its Kramers partner ($\bar{1}^-$).
As shown here, we put a bar on the orbital number to denote the Kramers partner.
Similarly, the size of the density matrix for the hydrogen molecule in this basis set is 8 ($N_D=4$).

We next explain the initial condition for the density matrix.
We choose the initial Heisenberg ket $|\Phi \rangle$ to be the ground states of the hydrogen atom and molecule,
which are obtained by the above explained computation method.
Namely, expressing $|\Phi \rangle = |\Phi_e \rangle \otimes  |\Phi_{ph} \rangle$,
where $|\Phi_e \rangle$ is the electron part and $|\Phi_{ph} \rangle$ is the photon part,
we use $|\Phi_e \rangle = \hat{e}^\dagger_{1^+}|0\rangle$ for the hydrogen atom,
and $|\Phi_e \rangle = \hat{e}^\dagger_{\bar{1}^+} \hat{e}^\dagger_{1^+}|0\rangle$ for the hydrogen molecule.
In general, the ground state of a $N_e$-electron system in the Hartree-Fock method is expressed as
$|\Phi_e \rangle = \prod_{i=1}^{N_e/2} \hat{e}^\dagger_{\bar{i}^+} \hat{e}^\dagger_{i^+}  |0\rangle$
when $N_e$ is even, and
$|\Phi_e \rangle = \hat{e}^\dagger_{((N_e+1)/2)^+} \prod_{i=1}^{(N_e-1)/2} \hat{e}^\dagger_{\bar{i}^+} \hat{e}^\dagger_{i^+}|0\rangle$ when $N_e$ is odd.
For later use, we here introduce the terminology ``occupied" orbitals.
If $|\Phi_e \rangle$ contains $\hat{e}^\dagger_{i^+}$, the $i$-th electron orbital is called ``occupied"
($i$ can be with or without bar),
and we can write $|\Phi_e \rangle = \prod_{i={\rm occupied}} \hat{e}^\dagger_{i^+}  |0\rangle$.
Then, using the anti-commutation relation,
the initial condition for the density matrix is
\begin{eqnarray}
{\cal E}_{n^a m^b}(t=0) &=& \left\{ \begin{array}{ll}
\delta_{nm} & (a=b=+,\, n:{\rm occupied}) \\
\delta_{nm} & (a=b=-) \\
0 & ({\rm otherwise}) \\
\end{array} \right. .
\end{eqnarray}
In particular, the vacuum expectation value at $t=0$, which is needed to compute Eq.~\eqref{eq:phys_rhoe}, is
\begin{eqnarray}
{\cal E}^0_{n^a m^b}(t=0)  &=& \left\{ \begin{array}{ll}
\delta_{nm} & (a=b=-) \\
0 & ({\rm otherwise}) \\
\end{array} \right. .
\end{eqnarray}

Other numerical details are as follows.
We work in the atomic units so that $m_e = e = \hbar = 1$, and $c = 137.035 999 679$.
The 1\,a.u. of time corresponds to 2.419$\times 10^{-17}$\,s or 24.19\,as.
As for the positions of the atomic nuclei, we locate them at the origin in the case of the hydrogen atom,
and at $(x, y, z)=(0,0,\pm0.7)$ in the case of the hydrogen molecule.
For both hydrogen 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}[{\cal E},\dot{\cal E}]$, Eq.~\eqref{eq:def_IjT}, 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.

\subsection{Effect of photon coherent state} \label{sec:coherent}
In this section, we show the results when the first approximation method explained
in Sec.~\ref{sec:approx_rad} is adopted.
Namely, we solve the time differential equation \eqref{eq:ev_dcalEdt} using Eq.~\eqref{eq:calO_1}.
As mentioned in the end of Sec.~\ref{sec:calc_setup}, the third term of Eq.~\eqref{eq:calO_1},
expressing the retarded potential, is neglected.
As for the photon initial Heisenberg ket, when there is no radiation field, $|\Phi_{ph} \rangle = | 0 \rangle $,
since $\langle \hat{a}_{\vec{p}_\sigma} \rangle = 0$, the forth term of Eq.~\eqref{eq:calO_1} is dropped.
In fact, $\langle \hat{a}_{\vec{p}_\sigma} \rangle = 0$ holds when $|\Phi_{ph} \rangle$ is any photon number state.
We may have $\langle \hat{a}_{\vec{p}_\sigma} \rangle \neq 0$ when $|\Phi_{ph} \rangle$ 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.~\cite{Loudon}, we denote it as $| \{ \alpha \} \rangle$.
This is characterized by the eigenvalues of the photon annihilation operator of each mode $(\vec{p},\sigma)$
as
\begin{eqnarray}
\hat{a}_{\vec{p}_\sigma} | \{ \alpha \} \rangle = \alpha(\vec{p},\sigma) | \{ \alpha \} \rangle,
\end{eqnarray}
where $\alpha(\vec{p},\sigma)$ is called spectral amplitude \cite{Loudon}.
In this paper, we use the delta-function type spectral amplitude and its center is chosen
to be a mode $(\vec{p}_j, \sigma_j)$, as
\begin{eqnarray}
 \alpha(\vec{p},\sigma) = \alpha_j \delta^{(3)}(\vec{p} - \vec{p}_j) \delta_{\sigma \sigma_j}. \label{eq:sa_deltafunc}
\end{eqnarray}
Such photon state corresponds to a classical oscillating electromagnetic field whose
propagating direction is $\vec{p}_j$, direction of circularly polarization is $\sigma_j$, and
amplitude is proportional to $\alpha_j$.
The period of the oscillations is determined by $p^0 = |\vec{p}_j|  $ as $2\pi/(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.~\ref{fig:H_coh}
and the case of the hydrogen molecule is in the upper panel of Fig.~\ref{fig:H2_coh}.
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 \times 10^{-4}$\,a.u.
As has been argued in Ref.~\cite{Ichikawa2013}, since this is very close to the period which
is determined from twice the mass of electron, $2\pi/(2m_e c^2)=1.67 \times 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" \cite{Ichikawa2013}.
In Ref.~\cite{Ichikawa2013}, we show this using the hydrogen atom.
In this paper, we show that the electron-positron 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.~\eqref{eq:sa_deltafunc}
with $\vec{p}_j/ |\vec{p}_j| =  (1,0,0)$ and $\sigma_j = +1$.
We compute the cases with $p^0 =  10$ and 20 for each hydrogen atom and molecule.
The case of $p^0=10$ (20) is shown in the middle (lower) panel in Figs.~\ref{fig:H_coh} and \ref{fig:H2_coh}.
In these panels, we can see that the electron-positron oscillations with the short period are
modulated by the longer period oscillations which are caused by the external oscillating
electromagnetic field.
In fact, the periods of the modulating oscillations seen in the panels for $p^0=10$ in both hydrogen atom and molecule
are close to $2\pi/(p^0 c)=4.59 \times 10^{-3}$\,a.u.
Similarly, in the panels for $p^0=20$, we see the periods of the modulating oscillations are close to
$2\pi/(p^0 c)=2.29 \times 10^{-3}$\,a.u.
We note that we have tuned the value of $\alpha_j$ for each case, in order to make
these effects visible clearly.
We have chosen $\alpha_j = 10^3$ ($10^4$) when $p^0=10$ (20) for the hydrogen atom
and $\alpha_j = 2 \times 10^4$ ($10^5$) when $p^0=10$ (20) for the hydrogen molecule.

We note that whether this ``electron-positron oscillations" is a real physical phenomenon
or an artifact of our model is an open question.
We can say with certainty that the electron-positron oscillations are caused by including the positron solutions
in our expansion functions for the Dirac field operator
(so they do not take place if we model the electron by the Schr\"{o}dinger field),
and, in fact, $(1^+,1^-)$-component of the density matrix oscillates in the case of the hydrogen atom simulation.
However, it is also true that the results are obtained with several approximations and
some ingredients of QED are missing.
In particular, including the retarded potential may affect the electron-positron oscillations.

\subsection{Effect of electron self-energy} \label{sec:selfenergy}
In this section, we show the results when the second approximation method explained
in Sec.~\ref{sec:approx_rad} is adopted.
Namely, we solve time differential equations \eqref{eq:ev_dcalEdt} and \eqref{eq:ev_deaedt}
using Eqs.~\eqref{eq:calO_2}, \eqref{eq:calP} and \eqref{eq:calQ}.
As in Sec.~\ref{sec:coherent}, 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, $|\Phi_{ph} \rangle = | 0 \rangle $,
in this section.
Therefore, the occupation number is zero for every photon mode $(\vec{p},\sigma)$, $n_{\vec{p}_\sigma}=0$.
However, even in such a case, as is expressed by the factor $(n_{\vec{p}_\sigma}+1)$ in the third term of Eq.~\eqref{eq:calP},
every photon mode contributes to the radiation term.
This is reasonable because,
as is mentioned in the end of Sec.~\ref{sec:qnum_excitation}, 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 $\vec{p}_\sigma$ in ${\cal E}_{n^a \vec{p}_\sigma m^b}$ to perform
numerical calculation.
We adopt the spherical coordinate system $(p^0, \theta, \phi)$ to express $\vec{p}$ and
use equally spaced grid points for each coordinate whose numbers are denoted by $N_{p^0}$, $N_\theta$ and $N_\phi$ respectively.
One more parameter we need to specify is the maximum of $p^0$, denoted by $p^0_{\rm max}$.
We first set $N_\theta=5$ and $N_\phi=4$, and compare the case with $(p^0_{\rm max}, N_{p^0}) = (10, 10)$ and
the case with $(p^0_{\rm max}, N_{p^0}) = (20, 20)$.
Since the results do not change, we adopt $(p^0_{\rm max}, N_{p^0}) = (10, 10)$ in the following.
As for the choice of $(N_\theta, N_\phi)$, the result with $(N_\theta, N_\phi) =(7, 8)$ is slightly different from
the case with $(N_\theta, N_\phi) =(5, 4)$, but it is almost same as the result with $(N_\theta, N_\phi) =(11, 10)$.
Therefore, in summary, we adopt $(p^0_{\rm max}, N_{p^0}, N_\theta, N_\phi) =(10, 10, 11, 10)$
as the photon-mode discretization parameters.

The results are shown in Fig.~\ref{fig:se}, the upper panel for the hydrogen atom and the lower panel for the hydrogen molecule.
As in Sec.~\ref{sec:coherent}, 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.~\ref{sec:coherent}.
However, in Fig.~\ref{fig:se}, 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 self-energy of the electron is the interaction energy between the electron and the electromagnetic field
 which is originated from the electron itself \cite{Jackson,SakuraiQFT}.
(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 $\hat{e}^\dagger \hat{a} \hat{e}$-type operator.
We would have a greater self-energy effect by considering the time evolution equation of
higher order operators such as the $\hat{e}^\dagger  \hat{a} \hat{a}^\dagger \hat{e}$-type operator,
which appears in the last term of Eq.~\eqref{eq:op_P}.
Another reason is that we have only included localized wavepackets for the expansion functions in Eq.~\eqref{eq:psi_expand2}.
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.

\section{Conclusion}  \label{sec: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 $\hat{e}^\dagger \hat{a} \hat{e}$-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 $\hat{e}^\dagger \hat{a} \hat{e}$-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.

\noindent
\section*{Acknowledgment}
This work is supported by Grant-in-Aid for Scientific research (No.~25410012)
from the Ministry of Education, Culture, Sports, Science and Technology, Japan.

redpajama_set_name: RedPajamaArXiv


Damascus steel isn't just for men!
One of our clients wanted a gift for his beloved to match his blackened Damascus steel wedding ring. We made a heart pendant for him to surprise his bride-to-be when he received his ring. This is a charming way to have your jewelry matching even if she doesn't have a Damascus ring.
Check our our reviews on WeddingWire!


"Professional. Fast. Explains well. Thanks a lot." -Brent and Jennifer.
My Family Obituary said: "Always available for help!"
He is a great teacher!
Gambler First Nation said: "He is a great teacher!"
Angusville General Store said: "Professional, Patient, Efficient, Met my Expectations...absoloutely excellent to work with!"
Richard took my ideas and created a website that was perfect for me!


"Absorbed in the ecstasy of the gopis, Lord Caitanya Mahaprabhu wished to see Lord Jagannatha in His original form as Krsna, the son of Nanda Maharaja, standing in Vrndavana and appearing very beautiful, His body curved in three places"
Absorbed in the ecstasy of the gopīs, Lord Caitanya Mahāprabhu wished to see Lord Jagannātha in His original form as Kṛṣṇa, the son of Nanda Mahārāja, standing in Vṛndāvana and appearing very beautiful, His body curved in three places. His desire to see that form was always increasing.
Facts about "Absorbed in the ecstasy of the gopis, Lord Caitanya Mahaprabhu wished to see Lord Jagannatha in His original form as Krsna, the son of Nanda Maharaja, standing in Vrndavana and appearing very beautiful, His body curved in three places"
This page was last modified on 15 September 2017, at 16:47.


If you are an experienced Welder, Structural Fitter, CNC Operator, Steel Detailer, Steel Erector, Field Installer, Labourer, Shipper Receiver, or Driver, then we would like to hear from you! We are willing to train individuals who are passionate, dedicated, and willing to learn.
BUILD your Career at Warwick Structures Group Ltd. Why should you consider us?
Investing in OUR future – with the construction of a new 212,000 square foot manufacturing facility – our new home in September 2015.


Cnbc Million Dollar Portfolio Challenge contestant beats 254,000 other stock traders and wins $500,000! Now you can literally copy his exact trades and safely make money month after month. Easily double your money! Stunning 92.4% accuracy.
Conversions Are Insane! Lethally Accurate Forex System Inc. System & Indicator Enabling Anybody To Make Shocking Amounts Of Money Trading The Forex Market With Only A Few Minutes A Day! Huge 70% Commission!
Membership Website For Traders. Converts Like Crazy.
Brand New! Killer Sales Page - High Converting Stocks Investment/make Money Site. Converts 1 : 20! Sales Are Flying In At The Higher Price! Top Affiliate Making Over $2000 Per Day! No Refunds!
Power3Marketing.com - 5.3% Conversion Rate Shows People Love Us!
Guaranteed revenue! Leverage the world's largest and fastest growing niche market! Our economy stinks but the medical industry and the elderly population are Booming! Invest in this 305 page ebook and discover everything to build a Nemt business.


text:
Our new series entitled 'From Zinc to RHEINZINK' begins with a small introduction to the all-rounder, zinc, as this trace element is not just good for our health. Zinc surrounds us everywhere in everyday life, and has been used as a building material for more than two centuries. Here, you can find out which processing options are available, how energy-efficient and sustainable the material is, and so much more.
All-rounder zinc
As a trace element, zinc is essential for human health; growth and wound-healing are positively influenced, too. Above, all, zinc is well-known for its role in our immune system, and the corresponding preparations are very popular, especially during the cold season. Moreover, zinc surrounds us everywhere in everyday life; from pressure-casting in fittings, and components in automobile construction, to hot-dip galvanised railings, to guttering, façade cladding or roofing, such as is manufactured at RHEINZINK.
Part 1: The material
Zinc has been used in construction for more than 200 years – titanium zinc is the modern version of this naturally-occurring metal. RHEINZINK is an alloy which consists of fine zinc and minimal quantities of titanium and copper. The fine zinc and titanium are sourced from a variety of countries of origin, such as Scandinavia, Spain, Belgium and Germany – rarely though from Mexico or Peru. Copper, as a recycling material, comes in the form of broken scrap for processing. The alloy additives are important for the patina colour and the special material-technological properties.
The alloy composition gives the materials not just the necessary strength, but the products can easily be shaped, too. These processing possibilities allow complicated details to be mastered and the realisation of any roof shapes. RHEINZINK titanium zinc, whether as cladding for a dormer, gable, balcony or chimney, represents an upgrade for the entire building due to its appealing and sophisticated appearance. Roof edges, porches and wall coverings shine in new splendour thanks to this noble material.
Nearly 100% of the titanium zinc used can be recycled with very little energy use and without any environmentally-harmful residues. The protective patina, which forms on the zinc surface over time, ensures a tremendously long life expectancy. Numerous historical buildings show that zinc roofs are completely intact, even after 200 years, and this without any notable maintenance effort. Furthermore, every RHEINZINK product already consists of 30% secondary material.
Part 2: Pre-alloying
The manufacturing process of RHEINZINK titanium zinc consists of several steps, which are carried out in the production at the Datteln plant. We take a look at the pre-alloying areas, the melting-casting-rolling line, the pre-weathering plant and then end with further processing. The first step is the pre-alloying.
Pre-alloying
The pre-alloying product area may only be accessed by skilled workers with precise instructions and special protective equipment as, in the truest sense of the word, it gets hot here! At a temperature of nearly 800°C, so-called pre-alloy ingots are produced from fine zinc, titanium and copper in a crucible induction furnace (capacity 2t). These ingots contain many times the titanium and copper content of the later rolled alloy.
It's getting hot!
What exactly is pre-alloying? This refers to alloys which are not used to manufacture end products, but are added to molten metal. The pre-alloy supports the process of metal bonding. For example, the variety of melting points make it very difficult to join various metals into a final rolled alloy. Thus, all three components of the alloy, as required by RHEINZINK, have different melting points: zinc 418°C; copper 1,083°C; titanium: 1,670°C.
The pre-alloy ingots are later added, in a certain ratio with fine zinc, to the large melting furnace, where a temperature of 'only' approx. 550°C is then required for melting.
Part 3: Melting-Casting-Rolling
The melting-casting-rolling line
In this chapter of our series entitled 'From Zinc to RHEINZINK', we present a really special feature. The RHEINZINK plant in Datteln is home to the only continuous melting-casting-rolling line worldwide!
In the 1960s, the RHEINZINK engineers played a vital role in the development of this patented melting-casting-rolling line, thus replacing the traditional packet rolling process. The difference to other manufacturers is that the special RHEINZINK alloy is processed in a continuous production process, from melting, to rolling, to coils of different metal thicknesses.
The first step is the melting. The ovens reach a temperature of 550°C and melt the pre-alloy together with the fine zinc ingots and/or return scrap into the final RHEINZINK alloy. The liquid alloy obtained is now introduced into the casting machine, one of the centrepieces of the RHEINZINK process. The technology upon which this is based creates the conditions for a uniform and, above, all, failure-free material quality. The circulating cooling water ensures that the RHEINZINK alloy is cooled so far below the melting point that a solid cast strand is produced.
The RHEINZINK method
Thereafter, five roll stands arranged in a row wait for the cast strand. In each roll stand, the thickness of the 15 mm thick strand is reduced by 35 to 50 percent. The final material thickness is between 0.5 - 2 mm. After the last roll stand, the now finished RHEINZINK coils are rolled up to become 20-tonne so-called large coils. The change takes place 'on the fly', as the continuous production process does not allow for any stoppages. This method is also a RHEINZINK development and unique worldwide. The still hot large coils (approx. 140°C) are then temporarily stored for at least two days to cool down to room temperature.
Part 4: The pre-weathering system
RHEINZINK coil is manufactured in a continuous process and, at the end, giant coils are created which gleam slightly and have a silvery appearance. We are referring to the bright rolled RHEINZINK material, as well as to our RHEINZINK CLASSIC product range. Over time, this surface forms a matt grey protective layer, the so-called patina. As the process of patina forming is more irregular, it sometimes has an uneven appearance, i.e. external factors, such as rain, influence the speed of the patina formation. Thus, bright rolled RHEINZINK, which is installed in a corner of a building protected against the rain, will naturally develop its patina layer slower than a roof-surface on the same building. It is somewhat comparable with wood. The table-top of a garden table usually weathers faster than the table-legs – especially the sides which slope inwards. The delayed patina formation is not a disadvantage, it is a natural process, but many architects, builders and developers prefer a uniform appearance right from the start. This is where our pre-weathering comes into play.
Since the 1980s, we have been offering pre-weathered surfaces. What we are speaking of here is our current prePATINA product range in the versions, blue-grey and graphite-grey. With these variants, the colour tone of the protective layer, i.e. the patina, is anticipated in our pre-weathering system. The natural surface characteristic of forming the life-long protective patina under environmental conditions is thereby retained. This process is also unique worldwide. The 'pre-weathered' products of our market competitors are either coated or painted/varnished and therefore no longer natural. The positive effect of the patina formation is therefore no longer present.
In our pre-weathering system, the bright rolled coils are subjected to pickling/staining and rinsed after a cleaning process. The whole procedure takes place in a closed, continuously-operating pickling/staining system, which meets the highest environmental requirements. Many measures are also taken that are far stricter than the current official environmental protection regulations: these include a dialysis system which converts the residual pickle/staining material into a reusable acid, thus drastically limiting the total acid requirement. In addition, a state-of-the-art catalyst system is used to clean the exhaust air, and the zinc sludge produced is recycled, too.
Whether the material becomes blue-grey or graphite-grey depends on the original alloy. For the graphite-grey material, the copper content is higher, for example. Incidentally, both bright rolled and pre-weathered materials have been fitted on our plant premises. Clearly, the Libeskind Villa was clad with prePATINA blue-grey. The façade elements at the zinc port have also been clad with blue-grey or granite-grey elements. In contrast, the zinc surfaces at the main entrance of the administration building were clad with bright rolled RHEINZINK in the 1960s. In the next part, we will take a closer look at the further processing department.
Part 5: Production // Further processing
Large coils have an impressive weight of about 20 tonnes and they are therefore naturally not viable for further processing. In the slitting/cut-to-length line, the coils are not only divided into smaller coils and sheets, but also subjected to the strength-bend-straightening process. This process pulls stresses out from the material which are created by the rolling process. The result is a wave-reduced material with straighter edges.
Up to this production step, the RHEINZINK semi-finished products are manufactured at the Datteln plant, i.e. coils and sheets, which can be further processed into finished products by the customers themselves. Naturally, RHEINZINK also has a wide product range, for example, it includes everything to do with roof drainage, or panels, and panels for roof and façade cladding.
This is why a large proportion of the clean-cut coils and sheets go on to further processing. There, automated processes meet manual work steps. Gutters, for example, are formed fully automatically, bundled and then packed onto pallets. The pipe-bending machine impressively bends a long straight pipe into a snake. After cutting, this bending creates individual pipe bends, but there are also areas where the employees reach for the soldering iron and connect components for roof drainage. The further processing department manufactures a total of more than 500 products in a wide variety of thicknesses and surface combinations. Some products are manufactured from our materials off-site by external, specialised suppliers, for example, RHEINZINK water-catch boxes are produced in southern Germany.
Small coils are also produced in the further processing department and prepared for transportation. Small coils weigh approx. 500 kg and are thus just 2.5 percent of the weight of the original large coil. After completion and packaging, all the products – whether gutters, pipes, panels, pipe clamps or small coils – are transferred to the logistics warehouse from where they are delivered to the customers.
Our little 'tour of the plant' ends with a visit to the further processing department. In the next part of our series, we will show you how our location in Datteln has developed since 1967.
Part 6: RHEINZINK history
Founded in 1967, the Datteln RHEINZINK-plant has evolved enormously in all areas over the decades.
Some of the milestones, apart from the construction of the rolling mill, are certainly the first steps in further processing in the early 1970s, the commissioning of the first pre-weathering plant at the end of the 1980s, and the successive further development of the office space, warehouse and production halls at the end of the 1990s and in the 2000s.
In 2006, the pre-weathering hall was opened in its current form. In 2009, the Libeskind Villa was built. Currently, a new large coil warehouse rounds off our factory premises in Datteln.
Have a look at our photos to get an idea of more than five decades of RHEINZINK history.
Part 7: The importance of roof drainage
There are numerous different requirements for both new buildings and renovations. Such requirements include the installation of a roof drainage system to safely drain rainwater and protect the building.
A roof drainage system is also indispensable for buildings as it protects the building's substance. RHEINZINK is a specialist for gutters, downpipes and a variety of accessories made of titanium zinc. The products are available in various finishes as well as in semi-circular or box-shaped designs. They therefore offer both a functional as well as a visually appealing solution to many different building styles.
When you take a closer look at the gutters of houses, the embossing reveals their origin. Gutters, downpipes or even end bases are embossed with a small RHEINZINK logo. But there are also roof drainage products from RHEINZINK without the tell-tale embossing. For example, reworkers use RHEINZINK coil material to manufacture their own gutters and pipes.
A quick overview of the advantages of a RHEINZINK-roof drainage system: precisely fitting system solutions for safe roof drainage, simple and cost-efficient installation, RHEINZINK - QUALITY ZINC and a 40-year warranty on our virtually maintenance-free material.
Part 8: The architectural details
Architectural details are elements that enrich the architectural impression of a building. A consistent overall appearance can ultimately only be achieved if these harmonize with façades or roof elements.
Complement with style
As an ideal complement, the RHEINZINK product range therefore includes a variety of façade details and architectural elements, which planners or architects can use to stylishly incorporate certain parts of a building into an overall concept, such as dormer and gable cladding, copings, roof edges or chimney head cladding.
Reliable and unique
There is a product for every architectural challenge in the RHIENZINK range. RHEINZINK flashings and claddings offer reliable protection from the elements on critical building parts, while at the same time impressing with their durability and giving the building that special something extra. What's more, they can also be customized to meet individual requirements.
Part 9: Façade claddings
A façade is an architectural design element; it defines the appearance of a building. However, it not only creates a representative appearance, but also has other important functions. For example, it protects the building from the weather, environmental influences and noise. RHEINZINK can also be used to masterfully set the scene for entire façades or partial surfaces.
Different installation systems
There are different installation systems, providing a very unique look for the building envelope. If you prefer a classic, handcrafted look, the angled standing seam with the typical optical interlocking joint or tiles in large or small formats are the perfect choice. If you prefer a more clean and straightforward look, you will certainly like the panel systems better. They can also be installed to resemble wood panels.
Appealing surfaces
Products are available in the various finishes depending on the respective system. Façade claddings made of RHEINZINK can be found on administrative buildings, public facilities, as well as private homes. The advantages: low-maintenance, easy to care for, long-lasting and simply a visual highlight.
Part 10: Roof coverings
It doesn't always have to be the classic roof tile; metal roofing is an attractive alternative, as it is lightweight, leak-proof and non-combustible. RHEINZINK has solutions for roofs that are weather-resistant, sustainable, long-lasting and have a highly appealing look. To achieve this, various installation techniques have been developed, which can be realized with RHEINZINK: the sheet and coil material is manufactured in further processing steps into either long panels or tiles in small or large formats.
Advantage of titanium zinc
The advantage of titanium zinc: Thanks to the handcrafted processing of the material, many shapes, curves and roof geometries are possible. Therefore, RHEINZINK can also be used for the entire building envelope. In addition, thanks to the various RHEINZINK surfaces - CLASSIC, prePATINA, artCOLOR, GRANUM and PRISMO, the right roofing material can be found for every building style and taste.
Part 11: RHEINZINK product lines
One of our core competencies is the development of innovative, high-quality products. The extensive knowledge and creativity of our development team, the latest trends in technology and design, and the needs of our customers are factors that ultimately lead from an idea to a product.
Some attractive examples from the extensive RHEINZINK product world are our premium surfaces: prePATINA, CLASSIC, artCOLOR, GRANUM and PRISMO, which are available in many designs and are suitable for a wide variety of purposes.
Durable and effective
RHEINZINK-CLASSIC and RHEINZINK-prePATINA are high-quality surfaces that over time form a blue-grey or graphite-grey patina. It is equipped with a natural self-regenerative effect and guarantees a long, maintenance-free life. RHEINZINK-artCOLOR is characterized by a high-quality PVDF coating and is also UV-resistant. The range encompasses seven effective colour shades, and colours can also be produced individually on request.
Noble and natural
GRANUM stands out with its matt and noble appearance. The phosphate coating is durable, weather-resistant and reveals the typical grain of the naturalpatina after installation. Even with the coloured glazed surface PRISMO, the typical zinc structure remains visible. A selection of six natural-looking colour shades is available.

redpajama_set_name: RedPajamaCommonCrawl


This leather loafer comes with a leather sole and oozes charisma. The upper has been embellished with a double-buckle detail and large kiltie flap in a matching tan colour. A classic versatile design defined by lean contours and energised by bold accents.
Leather loafer with a double buckle, large kiltie flap detailing and leather sole.


River herring is the common term used to describe two species of anadromous fish: alewives and blueback herring. Rivers along the eastern seaboard once teemed with these fish in the spring as they migrated from the ocean to breed in our rivers. Populations are in decline along the eastern seaboard. Humans' many dams and activities on land and in the sea are largely to blame. The federal government is considering listing these fish as endangered species.
River herring were recognized as being one of three "outstandingly remarkable values" in the Lamprey River's Wild and Scenic River designation by the US Congress. The continued presence of these fish in the Lamprey is important locally and nationally.
Click here to read the summary of a presentation by Kevin Sullivan of the New Hampshire Department of Fish and Game from December 8, 2015.
Click here for Kevin's slides.