We introduce a formalism to describe fractional charging of a molecule due to interactions with its environment. The interactions which induce fractional charging are contained in the Hamiltonian of the full system (molecule and environment). Such interactions can be singled out by expressing the Hamiltonian in a local spin orbital basis, and they are the main focus of this work. A reduced density operator for the molecule is derived starting from the Liouville-von Neumann equation for the full system by employing an explicitly constructed projection superoperator. By treating the molecule as an electronically open quantum system, we obtain a Redfield equation where the environment is included approximately. Phenomenological broadening of energy levels is included to mimic finite lifetimes of electronic states. The populations of the reduced density operator determine the mixture of different redox states and, hence, the fractional charging of the molecule. To illustrate the formalism, we use benzene physisorbed on a graphene sheet as a toy model. The work presented in this paper constitutes an initial step toward understanding molecules as electronically open quantum systems.

§ INTRODUCTION

State-of-the-art electronic structure models may describe a molecule interacting with its environment. Numerous approaches are developed to include such environmental effects, either quantum mechanically<cit.> or classically.<cit.> A common feature of these approaches is that they only include interactions that conserve the number of electrons in the molecule and environment separately. Hence, the molecule is treated as an electronically closed system.
However, if a full quantum-mechanical description of both the molecule and the environment is carried out, electrons are not strictly enforced to belong to either subsystem. Consequently, a full quantum-mechanical description allows for electron transfer between subsystems. The electron transfer may range from integer transfer (redox reaction) to partial transfer (fractional charging). Fractional charging is the manifestation of electrons being shared to some extent between the subsystems, and while integer electron transfer is widely studied<cit.>, fractional charging is largely ignored.

We will describe the fractional charging of a molecule due to the interaction with its environment. Through an analysis of the Hamiltonian of the full system in a basis of spin orbitals localized to either of the subsystems, the terms inducing fractional charging are identified. From these terms, an effective interaction operator is constructed. Our analysis of the Hamiltonian directly shows which one- and two-electron integrals are contained in the effective interaction strength.
We note that the effective operator is of the same form as in electron transport models<cit.>.
To describe the effect of the interaction on the molecule, we will use the density operator formalism and the theory of open quantum systems.<cit.> The theory of open quantum systems aims to describe the time evolution of a system in the presence of an environment, where only the environment's effect on the target system is included approximately. Previously, the theory has been used to model numerous phenomena, such as vibrational relaxation for spectroscopy<cit.> and energy transfer in electron transfer processes.<cit.>

The advantage of using density operators is that they can describe irreversible processes not attainable when using the electronic wave function. Examples are dissipation and decoherence<cit.>, and the extent of these effects depends on the interaction between the molecule and the environment.<cit.> Furthermore, the density operator allows for direct access to the state populations.
The time development of the reduced density operator describing only the molecule explicitly is derived from the Liouville-von Neumann equation for the full system. This is achieved by using the projection approach of Nakajima and Zwanzig<cit.>, resulting in two coupled equations. In the absence of the interaction operator, the two equations decouple and correspond to the Liouville-von Neumann equations for the separate subsystems. Hence, the interaction couples the electronic degrees of freedom in the molecule to the electronic degrees of freedom in the environment.

The coupled equations resulting from the Nakajima-Zwanzig approach are as challenging to solve as the Liouville-von Neumann equation for the full system. To simplify, the Born and Markov approximations are typically introduced.<cit.> These approximations assume that the environment is large and not significantly altered by the interaction, and that the environment's relaxation time is short. Additionally, assuming that the molecule-environment interaction is weak leads to the Redfield equation. Applying the secular approximation on top of this yields a master equation of Lindblad form, ensuring complete positivity of the dynamics<cit.>. The Redfield and Lindblad equations have been extensively used in quantum optics<cit.>, and for electron transport and charge separation in molecular system.<cit.>

We derive the Redfield equation using the exact eigenstates of the molecular electronic Hamiltonian for the non-interacting molecule and environment. Within a given one-electron basis, the exact eigenstates are the full configuration interaction (FCI) states.<cit.> It is important to note that due to the form of the interaction operator, we need to consider states with different number of electrons. Since FCI calculations are only applicable to very small molecular systems, we use Slater determinants as the many-electron basis for the reduced density operator in the derived form of the Redfield equation. The determinants do not provide a quantitative description of the molecule, but serve the purpose of illustrating the presented formalism. From propagating the equation to a steady state, the populations of the determinants are obtained.
The average number of electrons in the molecule may be fractional due to mixing of determinants with different number of electrons. In this formalism, fractional charging can be calculated as the deviation from the number of electrons of the isolated molecule.

In standard electronic structure theory, fractional orbital occupations have been extensively used to recover static correlation for isolated molecules.<cit.>. On the other hand, the particle-breaking Hartree-Fock (PBHF) model<cit.> is developed to describe the molecule as an electronically open quantum system. The PBHF wave function is obtained by using the time-independent Schrödinger equation and the variational principle. The interaction operator is parametrized prior to the optimization, and PBHF therefore represents a convenient approach to capture fractional charging effects. However, irreversible processes due to the interactions cannot be represented by the wave function.

The paper is organized as follows. In Section <ref>, we present the subsystem analysis of the Hamiltonian in the local spin-orbital basis and use FCI states of the non-interacting subsystems to construct the density operator for the full system. Next, we show the derivation of the Redfield equation for the populations of the reduced density operator. This involves the use of an explicitly constructed projection superoperator. In Section <ref>, we present an illustrative example by using a benzene molecule physisorbed on graphene as a toy model. We further discuss the developments necessary to obtain an appropriate electronic description of the molecule. Finally, Section <ref> provides a summary and concluding remarks.

§ THEORY

We start from the full configuration interaction (FCI) solution to the electronic problem of
a closed system, as this provides eigenstates of the Hamiltonian within the given one-electron basis. Next, we partition the closed system into two interacting subsystems, the molecule and its environment. The goal is to obtain a reduced description of the molecule while retaining the components that induce fractional charging. We construct the density operator for the full system (molecule and environment) starting from the exact solutions of the isolated subsystems. We then present the steps involved in the derivation of the Redfield equation, which gives the time evolution of the reduced density operator. Finally, we briefly discuss how to represent properties of the environment and the effect of energy level broadening.

§.§ Electronic structure theory for closed systems

The time-independent electronic Schrödinger equation is given by

\begin{align}
\hat H\Psi = E\Psi
\end{align}

where $\hat H$ is the molecular electronic Hamiltonian, $\Psi$ is the electronic wave function and $E$ is the energy.
The Hamiltonian in second quantization, excluding nuclear repulsion, is given by

\begin{align}
\label{eq:electronic_hamiltonian}
\hat H = \sum_{PQ}h_{PQ}\cre{P}\ann{Q} + \half\sum_{PQRS}g_{PQRS}\cre{P}\cre{R}\ann{S}\ann{Q}
\end{align}

where $\cre{P}$ and $\ann{P}$ are the creation and annihilation operators for spin-orbital $\phi_P$, respectively. The set of creation and annihilation operators for the spin orbitals obey the anticommutation rules $\comm{\cre{P}}{\cre{Q}}_+=\comm{\ann{P}}{\ann{Q}}_+=0$, and $\comm{\cre{P}}{\ann{Q}}_+=\delta_{PQ}$. Further, we have

\begin{align}
h_{PQ} &= \int \phi_P^*(\mathbf{x}_1)\Bigl(-\half\nabla^2-\sum_I\frac{Z_I}{|\mathbf{r}_1-\mathbf{R}_I|}\Bigr)\phi_Q(\mathbf{x}_1)d\mathbf{x}_1\\
g_{PQRS} &= \int \int \phi_P^*(\mathbf{x}_1)\phi_Q(\mathbf{x}_1)\frac{1}{|\mathbf{r}_1-\mathbf{r}_2|}\phi_R^*(\mathbf{x}_2)\phi_S(\mathbf{x}_2)d\mathbf{x}_1d\mathbf{x}_2\label{eq:two_el_integral}
\end{align}

The integration variable $\mathbf{x}_i$ includes spatial and spin coordinates of electron $i$,
$Z_I$ is the charge of nucleus $I$, and $\mathbf{r}_i$ and $\mathbf{R}_I$ are electron and nuclear coordinates, respectively. The Hamiltonian in eq <ref> represents a closed molecular system, and it commutes with the electron number operator $\hat N=\sum_P\cre{P}\ann{P}$. In other words, the number of electrons is conserved for the closed system.

The exact solution within a chosen atomic orbital basis is given by the FCI wave function,<cit.>

\begin{align}
\label{eq:fci_state}
\ket{\Psi_k} = \sum_{I}C_{Ik}\ket I
\end{align}

The sum over $I$ includes all determinants for distributing $N$ electrons among $M$ spin orbitals, with corresponding expansion coefficients $C_{Ik}$. The FCI states are orthogonal and may be normalized,

\begin{align}
\begin{split}
\braket{\Psi_k}{\Psi_l} = \delta_{kl}
\end{split}
\end{align}

The variational minimization of the energy with respect to the expansion coefficients amounts to diagonalizing $\mathbf{H}$, i.e., solving the FCI eigenproblem

\begin{align}
\mathbf{HC}_k=E_k\mathbf{C}_k
\end{align}

where $E_k$ is the energy of the $k$-th state, $\ket{\Psi_k}$, and the elements $H_{IJ}=\bra{I} \hat H \ket{J}$ are evaluated using the Slater-Condon rules.

§.§ Dividing the closed system into two interacting subsystems

We now consider an overall closed system, which we partition into a region of interest (the molecular system, $\systa$) and everything else (the environment, $\systb$). Henceforth, we will use the term environment for everything that is not treated explicitly. We further assume that we can transform the set of spin orbitals of the full system into a set where each spin orbital can be assigned to either the system or the environment. As long as there is no covalent bond between the subsystems, we know such localized spin orbitals exist.<cit.>
We note that this does not mean that the spin orbitals are all spatially localized, as they may be delocalized across either the molecule or the environment.

From now on, we will denote spin orbitals belonging to the molecule's orbital space with unbarred indices $\{\phi_p\}$, and orbitals belonging to the environment with barred indices $\{\phi_{\bar p}\}$, such that the set of all spin orbitals $\{\phi_P\} = \{\phi_p\}\cup\{\phi_{\bar p}\}$. We also use this distinction of indices for the second quantized operators. We may now separate the Hamiltonian in eq <ref> into three parts,

\begin{align}
\hat H = \hat H_\systa + \hat H_\systb + \hat H_{\systa\systb}
\end{align}

\begin{align}
\hat H_\systa &= \sum_{pq}h_{pq}\cre{p}\ann{q} + \half\sum_{pqrs}g_{pqrs}\cre{p}\cre{r}\ann{s}\ann{q}\\
\hat H_\systb &= \sum_{\bar p\bar q}h_{\bar p\bar q}\cre{\bar p}\ann{\bar q} + \half\sum_{\bar p\bar q\bar r\bar s}g_{\bar p\bar q\bar r\bar s}\cre{\bar p}\cre{\bar r}\ann{\bar s}\ann{\bar q}\\
\hat H_{\systa\systb} &= \hat H - \hat H_\systa - \hat H_\systb
\end{align}

$\hat H_\systa$ and $\hat H_\systb$ describe the two subsystems separately, and they commute. $\hat H_{\systa\systb}$ contains all interactions between the molecule and the environment. We may also divide the electron number operator into a number operator for the molecule and for the environment, $\hat N= \hat N_\systa + \hat N_\systb$. The subsystem number operators are

\begin{align}
\hat N_\systa &=\sum_p\cre{p}\ann{p}\\
\hat N_\systb &= \sum_{\bar p}\cre{\bar p}\ann{\bar p}
\end{align}

The molecule and environment number operators commute with their corresponding Hamiltonians, i.e., $\comm{\hat H_\systa}{\hat N_\systa} = \comm{\hat H_\systb}{\hat N_\systb} = 0$. In contrast, although the interaction Hamiltonian $\hat H_{\systa\systb}$ commutes with $\hat N$, only parts of it commute with $\hat N_\systa$ and $\hat N_\systb$ separately. We will therefore separate $\hat H_{\systa\systb}$ into two components,

\begin{align}
\hat H_{\systa\systb}=\hat H_{\systa\systb}^\text{pc}+\hat H_{\systa\systb}^\text{pb}
\end{align}

where $\hat H_{\systa\systb}^\text{pc}$ and $\hat H_{\systa\systb}^\text{pb}$ contain the particle-conserving and particle-breaking terms, respectively. By particle-breaking, we mean non-particle-conserving.<cit.>
These properties follow from the commutation relations

\begin{align}
\comm{\hat H_{\systa\systb}^\text{pc}}{\hat N_\systa} &= 0,\qquad
\comm{\hat H_{\systa\systb}^\text{pc}}{\hat N_\systb} = 0\\
\comm{\hat H_{\systa\systb}^\text{pb}}{\hat N_\systa} &\neq 0, \qquad
\comm{\hat H_{\systa\systb}^\text{pb}}{\hat N_\systb} \neq  0
\end{align}

Interactions contained in $\hat H^\text{pc}_{\systa\systb}$ may change the energies of the subsystems, but the respective particle numbers in each subsystem are unchanged. Numerous approaches<cit.> include these types of interactions and will therefore not be the focus of this work. However, interactions contained in $\hat H^\text{pb}_{\systa\systb}$ results in electron number fluctuations which may induce fractional charging. Here, we define fractional charging by

\begin{align}
\label{eq:fractional_charging}
\delta N_\systa = N_{\systa,0} - \langle \hat N_\systa \rangle
\end{align}

where $N_{\systa,0}$ is the number of electrons in the isolated system, and $\langle\hat N_\systa\rangle$ is the average number of electrons in the molecule upon interaction with the environment.

§.§.§ The particle-breaking interaction operator

We will now reduce $\hat H_{\systa\systb}^\text{pb}$ to an approximate effective interaction operator $\hat V$. All terms of $\hat H_{\systa\systb}^\text{pb}$ are

\begin{align}
\begin{split}
\hat H^\text{pb}_{\systa\systb}&= \sum_{p\bar q}\bigl(h_{p\bar q}\cre{p}\ann{\bar q} + h_{p\bar q}^*\cre{\bar q}\ann{p}\bigr)
+ \sum_{p\bar qrs}\bigl(g_{p\bar qrs}\cre{p}\cre{r}\ann{s}\ann{\bar q}
+  g_{p\bar qrs}^*\cre{\bar q}\cre{s}\ann{r}\ann{p} \bigr) \\
&+ \half\sum_{p\bar qr\bar s}\bigl(g_{p\bar qr\bar s}\cre{p}\cre{r}\ann{\bar s}\ann{\bar q} + g_{p\bar qr\bar s}^*\cre{\bar q}\cre{\bar s}\ann{r}\ann{p}\bigr) + \sum_{p\bar q\bar r\bar s}\bigl(g_{p\bar q\bar r\bar s}\cre{p}\cre{\bar r}\ann{\bar s}\ann{\bar q} + g_{p\bar q\bar r\bar s}\cre{\bar q}\cre{\bar s}\ann{\bar r}\ann{p}\bigr)
\end{split}
\end{align}

$\hat H^\text{pb}_{\systa\systb}$ will be reduced based on size of the charge distributions, $\phi^*_P(\mathbf{x}_i)\phi_Q(\mathbf{x}_i)$, which enter the two-electron integrals (see eq <ref>). Integrals such as $g_{p\bar q r \bar s}$ will generally be much smaller than integrals such as $g_{p\bar q  \bar r \bar s}$ and $g_{p\bar q r  s}$. This is due to the fact that $g_{p\bar q r \bar s}$ contains an integral over two charge distributions, $\phi^*_p(\mathbf{x}_1)\phi_{\bar q}(\mathbf{x}_1)$ and $\phi^*_r(\mathbf{x}_2)\phi_{\bar s}(\mathbf{x}_2)$, which both are necessarily small since $p$ and $r$ are localized to the system, whereas $\bar q$ and $\bar s$ are localized to the environment. In contrast, $g_{p\bar q r  s}$  contains an integral over the  charge distributions, $\phi^*_p(\mathbf{x}_1)\phi_{\bar q}(\mathbf{x}_1)$ and $\phi^*_r(\mathbf{x}_2)\phi_{s}(\mathbf{x}_2)$, of which $\phi^*_r(\mathbf{x}_2)\phi_{s}(\mathbf{x}_2)$ will be large, in particular for $r=s$.  Retaining only the dominant terms, and wrapping it into an effective one-electron operator, we obtain a form of $\hat V$ which connects Fock space states which differ in the number of electrons by one.

\begin{align}
\label{eq:interaction_approximated}
\hat V = \sum_{p\bar q}(V_{p\bar q}\cre{p}\ann{\bar q} + V_{p\bar q}^*\cre{\bar q}\ann{p})
\end{align}

The effective integral $V_{p\bar q}$ contains one-electron components as well as two-electron components of a charge distribution $\phi^*_p(\mathbf{x}_i)\phi_{\bar q}(\mathbf{x}_i)$ interacting with electrons in both the molecule and the environment.
$\hat V$ has the bilinear form often found in descriptions of bipartite open quantum systems<cit.>, but in this derivation, we have kept the physical content of it in terms of standard one- and two-electron integrals. However, note that although we know the content of  $V_{p\bar q}$, it can not be calculated without doing a full quantum mechanical calculation of the combined system and environment. It will therefore be parametrized based on knowledge of the magnitudes of the type of interactions it contains.

§.§ Density operator formalism for FCI states

We allow the two subsystems to interact by turning on the interaction $\hat V$ at some time $t=0$. We will now construct the density operator for the full system from the eigenstates of $\hat H_\systa$ and $\hat H_\systb$.
Using the effective interaction operator given in eq <ref>, the Hamiltonian for the interacting system $\systa+\systb$ is given by

\begin{align}
\hat H(t) = \hat H_\systa + \hat H_\systb + \hat V(t) \equiv  \hat H_0 + \hat V(t)
\end{align}

Turning on $\hat V(t)$ at time $t=0$ is enforced by using the Heaviside step function, $\hat V(t) = \theta(t)\hat V$. We regard $\hat H_0$ as the zeroth-order Hamiltonian and $\hat V(t)$ as a small perturbation. I.e., we assume that the interaction is weak, which is valid for the considered non-covalent interactions. We may construct common eigenstates for $\hat H_0$ from the eigenstates of $\hat H_\systa$ and $\hat H_\systb$ since $\hat H_0$ is additively separable.
The eigenstates for $\hat H_\systa$ and $\hat H_\systb$ are given by the FCI solutions

\begin{align}
\mathbf H_\systa\mathbf C^\systa_a &= E_a^\systa \mathbf C^\systa_a\\
\mathbf H_\systb\mathbf C^\systb_{\bar a} &= E_{\bar a}^\systb \mathbf C^\systb_{\bar a}
\end{align}

where the elements of the Hamiltonians are $\bra{I}\hat H_\systa\ket{J}$ and $\bra{\bar I}\hat H_\systb\ket{\bar J}$. The determinants are defined by $\ket{I}=\hat{A}_I^\dagger\vacket$ and $\ket{\bar I} = \hat{A}_{\bar I}^\dagger\vacket$, where $\hat{A}_I^\dagger$ and $\hat{A}_{\bar I}^\dagger$ are strings of creation operators for the system and environment, respectively.
We thus have states for the system, $\ket{a}$, and the environment, $\ket{\bar a}$,

\begin{align}
\ket{a} &= \sum_IC_{Ia}^\systa\ket{I}\label{eq:system_a_eigenstates}\\
\ket{\bar a} &= \sum_{\bar I}C_{\bar I\bar a}^\systb\ket{\bar I}\label{eq:system_b_eigenstates}
\end{align}

which satisfy $\braket{a}{b}=\delta_{ab}$ and $\braket{\bar a}{\bar b}=\delta_{\bar a\bar b}$. Using these solutions, we construct a combined basis for the system and environment

\begin{align}
\ket{a\bar a} = \sum_{I\bar I}C_{Ia}^\systa C_{\bar I\bar a}^\systb\hat{A}_{I}^\dagger\hat{A}_{\bar I}^\dagger\vacket
\end{align}

The orthonormality $\braket{a\bar a}{b\bar b}=\delta_{ab}\delta_{\bar a\bar b}$ follows from the properties of the FCI states (eqs <ref> - <ref>).
We note that the states $\ket{a\bar a}$ are not eigenstates of the full Hamiltonian $\hat H$, but they are eigenstates of $\hat H_0$,

\begin{align}
\hat H_0\ket{a\bar a} = \bigl(E_{a}^\systa + E_{\bar a}^\systb\bigr)\ket{a\bar a}
\end{align}

We build the density operator for the combined system from these zeroth-order solutions. In this basis, the density operator can be written as

\begin{align}
\rho(t) = \sum_{a\bar a b\bar b}\rho_{a\bar a,b\bar b}(t)\ket{a\bar a}\bra{b\bar b}
\end{align}

This form of the density operator will be used to construct a reduced density operator for the system, implicitly taking the environment into account.

§.§ The Liouville-von Neumann equation for the density operator

The time evolution of the density operator is governed by the Liouville-von Neumann equation, which in the interaction picture is given by

\begin{align}
\label{eq:lvn_interaction_picture}
\ddt\tilde\rho(t) = -i\comm{\tilde V(t)}{\tilde\rho(t)}
\end{align}

An operator $\hat \Omega$ in the interaction picture is given by the transformation $\tilde\Omega(t)=e^{i\hat H_0t}\hat\Omega e^{-i\hat H_0t}$. The goal is to obtain an equation for the time evolution of the molecular system $\systa$. To do this, we use the projection operator technique of Nakajima and Zwanzig<cit.>, where we define two orthogonal projection superoperators $\pproj$ and $\qproj$.
$\pproj$ is explicitly constructed such that the $\pproj$-projected density operator describes the dynamics of the molecule under the influence of the environment. This is achieved by

\begin{align}
\begin{split}
\label{eq:p_projection_operator}
\pproj\tilde\rho(t) &\equiv \sum_{a\bar a b}\ket{a\bar a}\bra{a\bar a}\tilde\rho(t)\ket{b\bar a}\bra{b\bar a}\\
&= \sum_{a\bar a b}\rho_{a\bar a,b\bar a}(t)e^{i\omega_{ab}t}\ket{a\bar a}\bra{b\bar a}
\end{split}
\end{align}

The exponential $e^{i\omega_{ab}t}$ results from the interaction picture, where $\omega_{ab}$ is the energy difference of the unperturbed molecule eigenstates, $\omega_{ab} = E_{a}^\systa-E_{b}^\systa$. From eq <ref>, we see that the $\pproj$ superoperator effectively averages out the environmental degrees of freedom. $\qproj$ is defined as the orthogonal complement of $\pproj$,

\begin{align}
\qproj\tilde\rho(t) &= (1-\pproj)\tilde\rho(t)
\end{align}

Since we only will consider weakly interacting subsystems, we approximate the elements $\rho_{a\bar a,b\bar a}(t)$ in eq <ref> as a product

\begin{align}
\label{eq:density_operator_decomposition}
\rho_{a\bar a,b\bar a}(t) \approx \sigma_{ab}(t)\rho_{\bar a \bar a}^\systb
\end{align}

where $\sigma_{ab}(t)$ are the elements of the reduced density matrix for the molecule. It is also assumed that the environment is large and has a short relaxation time relative to the system, such that the environment is unchanged by the interaction. Thus, the elements $\rho_{\bar a\bar a}^\systb$ are time-independent.
To get a reduced description of the molecule, a decomposition of the density operator $\rho(t)$ is essential.<cit.>

We can construct two coupled equations by projecting the Liouville-von Neumann equation (eq <ref>) with $\pproj$ and $\qproj$,

\begin{align}
\ddt\pproj\tilde\rho(t) &= -i\pproj\comm{\tilde V(t)}{\pproj\tilde\rho(t)} - i\pproj \comm{\tilde V(t)}{\qproj\tilde\rho(t)}\label{eq:p_projected_liouville}\\
\ddt\qproj\tilde\rho(t) &= -i\qproj\comm{\tilde V(t)}{\pproj\tilde\rho(t)} - i\qproj \comm{\tilde V(t)}{\qproj\tilde\rho(t)}\label{eq:q_projected_liouville}
\end{align}

where we also have inserted the identity $1=\pproj+\qproj$ in front of the $\tilde\rho(t)$ inside the commutator. Combined, these equations are still equivalent to the Liouville-von Neumann equation for $\tilde\rho(t)$ and are as challenging to solve as the full problem. However, they serve as a convenient starting point for constructing an equation for the elements $\sigma_{ab}(t)$ of the reduced density matrix for the molecule.

§.§ The Redfield equation for the reduced density matrix

We will now introduce approximations into eqs <ref> - <ref>. A formal integration of eq <ref> gives

\begin{align}
\begin{split}
\label{eq:q_projected_liouville_integrated}
\qproj\tilde\rho(t) = -i\int_0^t\qproj\comm{\tilde V(t')}{\pproj\tilde\rho(t')}dt'
-i\int_0^t\qproj\comm{\tilde V(t')}{\qproj\tilde\rho(t')}dt'
\end{split}
\end{align}

where $\qproj\tilde\rho(0)=0$ since the interaction $V(t)$ is turned on at time $t=0$. Following the steps in the derivation of the Redfield equation<cit.>, we assume the weak coupling limit, i.e., the interaction $\hat V$ is small compared to the energies of the molecule $E_a^\systa$. We approximate the above equation to the lowest order in $\hat V$ by disregarding the last integral, such that $\qproj\tilde\rho(t)$ no longer depends on itself. We further assume that the environment is large and unchanged by the interaction $\hat V$.
Inserting eq <ref> (excluding the omitted integral) back into eq <ref>, we obtain

\begin{align}
\label{eq:redfield_equation_double_commutator}
\ddt\pproj\tilde\rho(t) = -\int_0^t\pproj\comm{\tilde V(t)}{\comm{\tilde V(t')}{\pproj\tilde\rho(t')}}dt'
\end{align}

where we have used that $\pproj \tilde V(t) \pproj\tilde\rho(t)=0$. By further assuming that memory effects are negligible, we get in the Markovian limit the following equation for $\sigma_{ab}(t)$,

\begin{align}
\label{eq:redifeld_equation_general}
\ddt\sigma_{ab}(t) = -i\omega_{ab}\sigma_{ab}(t) +\sum_{cd}R_{abcd}\sigma_{cd}(t)
\end{align}

$\mathbf{R}$ is the Redfield tensor defined through the tensors $\boldsymbol\Gamma^+$ and $\boldsymbol\Gamma^-$,

\begin{align}
\begin{split}
R_{abcd} = &(-1)^{\Delta N_{cd}}\bigl(\Gamma_{dbac}^+(\omega_{da}) + \Gamma_{dbac}^-(\omega_{bc})\\
\delta_{bd}\sum_e\Gamma_{aeec}^+(\omega_{de})-\delta_{ac}\sum_e\Gamma_{deeb}^-(\omega_{ec})
\bigr)
\end{split}
\end{align}

$\Delta N_{cd}$ is the difference in the numbers of electrons between states $\ket{c}$ and $\ket{d}$. The $\boldsymbol\Gamma^x$ tensors ($x=+,-$) are defined by

\begin{align}
\begin{split}
\Gamma_{abcd}^x(\omega) = \sum_{pq}
&\bra{a}\cre{p}\ket{b}\bra{c}\ann{q}\ket{d}\int_0^\infty  e^{i\omega\tau}C_{pq}(x\tau)d\tau\\
\sum_{pq}&\bra{a}\ann{q}\ket{b}\bra{c}\cre{p}\ket{d}\int_0^\infty  e^{i\omega\tau}D_{pq}(x\tau)d\tau
\end{split}
\end{align}

By introducing the notation $\langle\cdot \rangle_\systb$ for the environment average (see Supporting Information for details), we define the time-correlation functions $C_{pq}(\tau)$ and $D_{pq}(\tau)$,

\begin{align}
C_{pq}(\tau) &= \sum_{\bar p\bar q}V_{p\bar p}V_{q\bar q}^*\bigl\langle\tann{\bar p}(\tau)\tcre{\bar q}(0)\bigr\rangle_\systb\label{eq:C_correlation_function}\\
D_{pq}(\tau) &= \sum_{\bar p\bar q}V_{p\bar p}V_{q\bar q}^*\bigl\langle\tcre{\bar q}(\tau)\tann{\bar p}(0)\bigr\rangle_\systb\label{eq:D_correlation_function}
\end{align}

Since we are interested in steady-state properties (long time scales) of the reduced density matrix, we continue by assuming the secular approximation in eq <ref>. This guarantees a master equation of Lindblad form, thereby overcoming any potential positivity issues of the Redfield equation.<cit.> The secular approximation exploits that off-diagonal elements of the reduced density matrix either oscillate rapidly compared to the evolution of the diagonal elements or decay to zero. Hence, the off-diagonal elements may be disregarded.<cit.>
Eq <ref> then simplifies to

\begin{align}
\label{eq:redifeld_equation_secular}
\ddt\sigma_{aa}(t) &= \sum_b k_{a\leftarrow b}\sigma_{bb}(t)
- \sum_b k_{b\leftarrow a}\sigma_{aa}(t)
\end{align}

where the transfer rates are defined by (see Supporting Information)

\begin{align}
\begin{split}
\label{eq:population_transfer_rate}
k_{b\leftarrow a}
= \sum_{pq}
&\bra{a}\cre{p}\ket{b}\bra{b}\ann{q}\ket{a}\int_{-\infty}^\infty  e^{i\omega_{ab}\tau}C_{pq}(\tau)d\tau\\
\sum_{pq}&\bra{a}\ann{q}\ket{b}\bra{b}\cre{p}\ket{a}\int_{-\infty}^\infty  e^{i\omega_{ab}\tau}D_{pq}(\tau)d\tau
\end{split}
\end{align}

Eq <ref> now has the form of a Pauli master equation.<cit.> By collecting the diagonal elements (populations) in a vector $\mathbf{P}$ with elements $P_a=\sigma_{aa}$, and the transfer rates in a matrix $\mathbf{W}$ with off-diagonal elements $W_{ab}=k_{a\leftarrow b}$ and diagonal elements $W_{aa}=-\sum_bk_{b\leftarrow a}$, we can write the equation as

\begin{align}
\label{eq:pauli_master_equation}
\ddt \mathbf{P} = \mathbf{WP}
\end{align}

Due to the form of the diagonal elements $W_{aa}$, the matrix $\mathbf{W}$ will always contain linearly dependent rows. This entails that at least one steady-state solution can be found by solving for its null space.

§.§ A mean-field approximation to the Redfield equation

The above derivation assumes that we have the FCI states for the non-interacting subsystems, and the introduced approximations so far have been those of the Redfield equation. We will now start by making approximations for the environment by assuming mean-field states and a continuous spectrum. The correlation functions in eqs <ref> - <ref> can be simplified to (see Supporting Information)

\begin{align}
\label{eq:correlation_function_minus_simplified}
C_{pq}(\tau) &= \int d\omega V_{p}(\omega)V_{q}^*(\omega)e^{-i\omega\tau}[1-\bar n(\omega)]\\
\label{eq:correlation_function_plus_simplified}
D_{pq}(\tau) &= \int d\omega V_p(\omega)V_q^*(\omega)e^{i\omega\tau}\bar n(\omega)
\end{align}

where $\bar n(\omega)$ is the average number of electrons at the given energy $\omega$, and $V_p(\omega)$ is the continuous representation of $V_{p\bar p}$. For a given finite temperature $T$ and a chemical potential $\mu$, the average number of electrons can, e.g., be given by the Fermi-Dirac distribution,

\begin{align}
\bar n(\omega) = \Bigl[\exp\Bigl(\frac{\omega-\mu}{k_BT}\Bigr)+1\Bigr]^{-1}
\end{align}

with $k_B$ being the Boltzmann constant. By inserting the simplified correlation functions in eqs <ref> - <ref> into the expression for the population transfer rates in eq <ref> we obtain

\begin{align}
\begin{split}
\label{eq:population_transfer_rate_2}
k_{b\leftarrow a}
= \sum_{pq}
\bra{a}\cre{p}\ket{b}\bra{b}\ann{q}\ket{a} J_{pq}(\omega_{ab})[1-\bar n(\omega_{ab})]\\
\bra{a}\ann{q}\ket{b}\bra{b}\cre{p}\ket{a} J_{pq}(\omega_{ba})\bar n(\omega_{ba})
\end{split}
\end{align}

where we have introduced $J_{pq}(\omega) \equiv 2\pi V_p(\omega)V_q^*(\omega)$, which we will refer to as the spectral density function.

We now approximate the electronic states of the molecule by using single determinants expressed in the same spin orbital basis. The transfer rates in eq <ref> can then be simplified by noting that $\bra{a}\cre{p}\ket{b}\bra{b}\ann{q}\ket{a}=0$ if $p\neq q$. Since the states $\ket{a}$ and $\ket{b}$ are expressed in the same spin orbital basis, the energy difference is $\omega_{ab}=\varepsilon_p$, where $\varepsilon_p$ is the orbital energy of $\phi_p$.
By writing $J_p(\omega) \equiv J_{pp}(\omega) = 2\pi|V_p(\omega)|^2$, the resulting expression for the population transfer rates is

\begin{align}
\label{eq:transfer_rate_single_determinants}
\begin{split}
k_{b\leftarrow a}
= \sum_{p}\Bigl(
\bra{a}\cre{p}\ket{b}\bra{b}\ann{p}\ket{a} J_{p}(\eps_p)[1-\bar n(\eps_p)]\\
\bra{a}\ann{p}\ket{b}\bra{b}\cre{p}\ket{a} J_{p}(\eps_p)\bar n(\eps_p)
\Bigr)
\end{split}
\end{align}

In principle, one may include all determinants. However, the choice of which determinants to include may be motivated by the fact that high-energy (relative to the chemical potential of the environment) spin orbitals will remain unoccupied. Similarly, low-energy spin orbitals will remain occupied. Spin-orbitals with energies close to the chemical potential therefore constitute an active space for constructing relevant determinants for the reduced density matrix.

§.§.§ The spectral density function

The spectral density function $J_p(\omega)$ contains information about the environment's density of states and the interaction strength of the system orbitals with the environment.
In this work, we approximate the interaction strengths to be equal for all spin orbitals. In this case, the spectral density function becomes proportional to the environment's electronic density of states. The density of states can either be modeled or measured experimentally.<cit.> We will estimate the density of states for the considered molecular environment by calculating its relative number of orbitals per some fixed energy spacing. In the wide band approximation, $J_p(\omega)$ is set to a constant. This is a reasonable choice for metallic environments such as gold or silver with a near-constant density of states around the Fermi energy.<cit.>

§.§.§ Energy level broadening

Excited states have finite lifetimes, and therefore broadened energy levels, due to interactions not captured by the electronic Hamiltonian.<cit.> In electronic structure theory, this has been accounted for by introducing phenomenological damping.<cit.> This relaxes the resonance requirement for transitions between molecule states.
We incorporate this effect by replacing Dirac delta functions in the time-correlation functions with normalized Gaussians with standard deviations $\gamma$ (see Supporting Information). Eq <ref> is therefore modified to

\begin{align}
\label{eq:transition_rate_with_energy_level_broadening}
\begin{split}
k_{b\leftarrow a}
= \frac{1}{\sqrt{2\pi\gamma^2}}\sum_{p}\Bigl(&\bra{a}\cre{p}\ket{b}\bra{b}\ann{p}\ket{a}\int_{-\infty}^\infty J_{p}(\omega)\exp\Bigl(-\frac{(\omega_{ab}-\omega)^2}{2\gamma^2}\Bigr)[1-\bar n(\omega)]d\omega\\
&\bra{a}\ann{p}\ket{b}\bra{b}\cre{p}\ket{a}\int_{-\infty}^\infty J_{p}(\omega)\exp\Bigl(-\frac{(\omega_{ba}-\omega)^2}{2\gamma^2}\Bigr)\bar n(\omega)d\omega
\Bigr)
\end{split}
\end{align}

We have now arrived at the final expression for the population transfer rates.

§ RESULTS AND DISCUSSION

In this section, we illustrate how the Redfield equation can be used to describe the fractional charging of a molecular system upon interaction with an environment. We note that we include only the novel particle-breaking interactions described in Section <ref>. However, particle-conserving interactions must also be considered for realistic chemical applications. As a model system, we consider a benzene molecule physisorbed on a finite graphene surface. The benzene molecule is the target for which we create a reduced density matrix, while the graphene sheet is the environment. We will explore how the fractional charging of benzene depends on the chemical potential, temperature, and energy level broadening.
In Section <ref>, we describe necessary computational details, while results are presented in Section
<ref>. Future prospects of using the proposed formalism are discussed in Section

§.§ Computational details

We will use a restricted Hartree-Fock calculation of isolated, neutral benzene (see Supporting Information for geometry specification) using the cc-pVTZ basis<cit.> as a reference determinant. All determinants will be constructed from these spin orbitals.
An active space of five molecular orbitals (i.e., ten spin orbitals) is chosen, see Figure <ref>. By distributing between zero and ten electrons among these orbitals, we obtain 32 possible determinants. The energy differences between the determinants are given by $\omega = \pm\varepsilon_p$, with $+$ for the addition and $-$ for the removal of an electron as mentioned in Section <ref>.

We describe the graphene surface, which plays the role of the environment, in a simplified manner. We assume that the graphene is unaffected by the interaction with the benzene molecule, as motivated by the many degrees of freedom (see Section <ref>). Further, all spin orbitals are assumed to interact with the same strength. The spectral density function is estimated by a Hartree-Fock calculation on a finite graphene sheet using the cc-pVDZ basis<cit.> and an energy spacing of 0.05 a.u. The density of states is calculated for sheets consisting of between 96 and 216 carbon atoms. Increasing the size of the graphene sheet does not change the density of states considerably, see Figure <ref>. We will therefore use the sheet containing 216 carbon atoms as the environment. The functional form resembles the reported density of states obtained from a density functional theory calculation presented in Ref. .

The Redfield equation (eq <ref>) using the population transfer rates of eq <ref> is propagated in time by Euler integration<cit.> until steady state is reached. The initial condition for the time propagation is the neutral benzene molecule.

Energy level diagram for the two degenerate highest occupied molecular orbitals and the three lowest unoccupied molecular orbitals of the neutral benzene molecule. The illustrated electron configuration corresponds to the ground state. The orbital energies are, in atomic units, $\eps_\text{HOMO} = -0.336252$, $\eps_\text{LUMO} = 0.129394$, and $\eps_{\text{LUMO}+1} = 0.145115$.

The approximated density of states calculated for graphene surfaces of different sizes. The density of states was calculated to be the relative number of orbitals within a window of 0.05 a.u. The orbitals were obtained from restricted Hartree-Fock calculations in the cc-pVDZ basis.

§.§ Numerical illustration: Fractional charging of benzene at steady state

Here, we show how the presented model can induce fractional charging of a molecule. By using the steady-state solution of benzene, we illustrate how fractional charging depends on the chemical potential of the environment and energy level broadening. We note that realistic temperatures will not impact fractional charging. If the energy differences between the molecule's charge states are on the order of $\SI{1}{eV}$, inducing fractional charging would require temperatures around $\SI{e4}{\kelvin}$. At such elevated temperatures, our molecule would break down. Therefore, we set the temperature to a realistic constant, $T=\SI{300}{\kelvin}$, in our calculations.

First, we use the wide band approximation for the spectral density function, $J_p(\omega)$. The results are shown in Figure <ref>. When the broadening $\gamma$ is small (close to zero) there is little to no fractional charging. In this case, the main effect is integer charging when the chemical potential crosses the energy differences of the integer charge states. When the chemical potential matches exactly the energy differences, the populations of the corresponding states are equal.
The values for $\delta N_\systa$ at $\mu=0$ are listed in Table <ref>, where the dependence of $\delta N_\systa$ on the broadening is seen. Even at $\mu=0$, significant fractional charging effects are seen for $\gamma > 0.03$. By significant, we mean fractional charging effects that influence the energy on the order of $10^{-3}$ Hartree.
We note that the fractional charging is zero when the chemical potential equals the midpoint between the HOMO and LUMO energy levels (see Figure <ref>), regardless of the broadening. The symmetry around this point is due to the symmetry of the wide band approximation.

Figure <ref> shows the fractional charging as a function of the chemical potential where the calculated electronic density of states is used as the spectral density function. Using the density of states shows a qualitative difference compared to using the wide band approximation. Notably, at $\mu=0$, the same energy level broadening gives a less significant fractional charging (see Table <ref>). This is because the density of states (see Figure <ref>) is close to zero around $\mu=0$. Further, the curves are no longer symmetric about neither the energy differences between the charge states nor about the midpoint between the HOMO and LUMO energies. These asymmetries are a direct consequence of the asymmetry of the density of states. The qualitative differences show that the choice of the spectral density function and the broadening and how it is modeled significantly affects the properties of the system of interest.

Broadening $\gamma$ [a.u.]     $\delta N_\systa$ for WBA     $\delta N_\systa$ for DOS

0.10	         0.537         0.167

0.09	         0.409         0.110

0.08	         0.283         0.064

0.07	         0.169         0.031

0.06	         0.080         0.001

0.05	         0.026         0.000

0.04	         0.003         0.000

0.03	         0.000         0.000

0.02	         0.000         0.000

0.01	         0.000         0.000

0.00	         0.000         0.000

Fractional charging $\delta N_\systa$ of benzene at zero chemical potential for different values of the energy level broadening $\gamma$. The results are shown for the two models for the spectral density function; the wide band approximation (WBA) and the calculated electronic density of state for graphene (DOS).

Fractional charging as a function of the chemical potential of the environment at different energy-level broadenings $\gamma$ with the wide band approximation as the spectral density function. The region marked in a) is enlarged in b). The temperature was set to $T=\SI{300}{\kelvin}$.

Fractional charging as a function of the chemical potential of the environment at different energy-level broadenings $\gamma$ with the calculated electronic density of states for graphene as the spectral density function. The region marked in a) is enlarged in b). The temperature was set to $T=\SI{300}{\kelvin}$.

§.§ Future prospects of using the Redfield equation for electronically open molecular systems

The presented results are based on several key approximations, which we will now address. These include the assumptions of the interaction operator, the system states, and the inherent assumptions of the Redfield equation. The approximations will be discussed in the context of describing a molecular system interacting weakly with an environment, and with the electronic structure of the molecule as the main objective. First, particle-conserving interactions (see Section <ref>) must also be included to appropriately account for the interaction with the environment. The particle-breaking interactions allow for describing fractional charging, but as discussed previously, we expect them to be only a small correction to the total interaction. The effect of particle-conserving interactions can, for example, be included by well-established quantum mechanical embedding methods or by a classical atomistic approach<cit.>. Second, we have used single determinants with a common set of spin orbitals as the basis for the system's reduced density matrix. Knowing that single-determinantal descriptions lack electron correlation, using correlated wave function models in the Redfield equation must be explored. A correlated model (e.g., the configuration interaction or coupled cluster wave functions) would change the Redfield equation's dynamics and result in a different steady-state solution. Third, the optimal spin orbitals change when we get a mixture of multiple charge states. However, relaxing the spin orbitals while propagating the Redfield equation will greatly increase the complexity and computational cost of the method. In addition, recomputing spin orbitals for the fractionally charged molecule would require a wave function method that handles non-integral electron numbers and fractional occupations. However, we expect the effect of fractional charging to be small such that the spin orbitals of the neutral molecule dominate.

The approximations for the Redfield equation may also be explored to assess its validity for describing fractional charging. The Born approximation is appropriate for non-covalent interactions between the system and a large environment that is unchanged by the interaction. Consequently, the Redfield equation is less suitable for describing environments consisting of finite-sized molecules. For smaller environments, memory effects become more significant. Further, if the system and the environment are similar in size, the Markov approximation breaks down due to the inability to separate the time scales.<cit.> Examining the validity of the Born and Markov approximations would require comparisons to more sophisticated methods, such as the hierarchical equations of motion.<cit.>

Another simplification is the assumption of equal interaction strengths for all benzene orbitals. In reality, interactions vary with distance and energy, and molecular orbitals are not uniquely defined. However, relative interaction strengths must be modeled in some way since explicit treatment of the environment is undesired. In principle, calculation of electronic coupling elements can be achieved by methods such as constrained DFT<cit.> or generalized Mulliken-Hush.<cit.>

§ SUMMARY AND CONCLUSIONS

In this work, we introduce a reduced density operator formalism to include the effects of fractional charging of a molecule due to a non-covalent interaction with its environment.  In order to develop the formalism, we start by analyzing the full configuration interaction solution in a spin orbital basis which is localized to either the molecule or the environment. In this basis, the terms of the molecular electronic Hamiltonian for the full system may be categorized into a molecule Hamiltonian, an environment Hamiltonian, and an interaction Hamiltonian. The interaction Hamiltonian may further be divided into particle-conserving and particle non-conserving (particle-breaking) terms by considering the commutators with the particle number operator for the molecule. The particle-breaking interactions and the fractional charging they induce are the focus of this work, and the paper thus addresses a type of interaction that has largely been neglected by the quantum chemistry community.

We include the particle-breaking interactions as a perturbation to the non-interacting molecule and environment. By introducing an explicit form of the projection superoperator  we may use the Nakajima-Zwanzig formalism from the theory of open quantum systems. This allows us to derive an equation for the reduced density operator of the molecule starting from the Liouville-von Neumann equation for the full density operator. We employ the Born and Markov approximations to obtain a Redfield equation for the populations of the reduced density operator,  where the effect of the environment is included approximately. The effect of the environment is modeled through the spectral density function. In this work, it is taken to be either a constant (the wide band approximation) or the computed electronic density of states for the environment. The fractional charging results from non-zero populations for molecular states with different integer charges and signifies that electrons are, to a certain extent, shared with the environment.  Further, a phenomenological broadening of energy levels is introduced to mimic the finite lifetimes of electronic states. We present numerical illustrations for a mean-field approach of the formalism, using benzene physisorbed on graphene as a toy model.

The work presented in this paper represents an important step toward understanding molecules as electronically open quantum systems. To achieve quantitative accuracy, it will in the future be necessary to include particle-conserving interactions and use correlated electronic states as the many-electron basis.

§ ACKNOWLEDGEMENTS

We thank Associate Professor Thorsten Hansen for fruitful discussions regarding Redfield theory.
I-M.H and B.S.S. acknowledge funding from the Research Council of Norway through
FRINATEK project 325574. I-M.H acknowledges support from the Centre for Advanced Study in Oslo, Norway, which funded and hosted her Young CAS Fellow research project during the academic year of 22/23 and 23/24. J.P. acknowledges financial support from the Technical University of Denmark within the Ph.D. Alliance Programme, Knud Højgaards Fond (grant no. 23-02-2571), William Demant Fonden (grant no. 23-3137), and Ingeniør Alexandre Haynman og fru Nina Haynmans Fond (grant no. KR-4110605).

get arXiv to do 4 passes: Label(s) may have changed. Rerun