Quantum Statistical Mechanics of Electronically Open Molecules: Reduced Density Operators

TL;DR

This paper introduces a new way to describe the quantum behavior of molecules interacting with their environment by accounting for electron sharing.

Contribution

A novel reduced density operator is proposed that includes particle-breaking interactions and resolves ambiguities in fermionic partial trace operations.

Findings

A new reduced density operator is defined for electronically open molecules.

The framework allows explicit electron occupancy in the environment at any theoretical level.

The approach is compatible with various environment models and includes a generalized chemical potential.

Abstract

We present a reduced density operator for electronically
open molecules
by explicitly averaging over the environmental degrees of freedom
of the composite Hamiltonian. Specifically, we include the particle-number
nonconserving (particle-breaking) interactions responsible for the
sharing of electrons between the molecule and the environment, which
are neglected in standard formulations of quantum statistical mechanics.
We propose an unambiguous definition of the partial trace operation
in the composite fermionic Fock space based on composite states in
a second quantization framework built from a common orthonormal set
of orbitals. Thereby, we resolve the fermionic partial trace ambiguity.
The common orbital basis is constructed by spatial localization of
the full orbital space, in which the full composite Hamiltonian naturally
partitions into a molecule Hamiltonian, an environment Hamiltonian,
and an interaction Hamiltonian. The new reduced density operator is
based on the approximation of commutativity between the subsystem
Hamiltonians (i.e., molecule and environment Hamiltonians) and the
interaction Hamiltonian, but our methodology provides a hierarchical
approach for improving this approximation. The reduced density operator
can be viewed as a generalization of the grand canonical density operator.
We are prompted to define the generalized chemical potential, which
aligns with the standard interpretation of the chemical potential,
apart from the possibility of fractional rather than strictly integer
electron transfer in our framework. In contrast to standard approaches,
our framework enables an explicit consideration of the electron occupancy
in the environment at any level of theory, irrespective of the model
used to describe the molecule. Specifically, our reduced density operator
is fully compatible with all possible level-of-theory treatments of
the environment. The approximations that render our reduced density
operator identical to the grand canonical density operator are (i)
restriction of excitations to occur within the same orbitals and (ii)
assumption of equal interaction with the environment for all molecule
spin orbitals (i.e., the wide-band approximation).

Figures

Illustration of our partitioning of a canonical molecular system energetically equilibrated with its surroundings into two interacting subsystems, namely, the molecule and the environment, which can exchange electrons.

Direction of electron flow in non-equilibrium situations can be predicted by the average environment orbital occupation. The different situations are favored electron transfer to the environment (left), equilibrium between electron addition and removal (center), and favored electron transfer to the molecule (right). The molecule’s ability to donate or accommodate an electron is taken into account by the molecule number operator. Therefore, the direction of electron flow is independent of the specific molecule spin orbitals.

Full text

Introduction

1

An open quantum system exchanges energy and/or particles with its environment, ?−? ? ? ? ? ? ? ? ? and we advocate viewing molecules as open quantum systems. Meanwhile, we are not interested in explicitly describing the composite system (molecule and environment). Instead, we pursue reduced descriptions aiming at explicitly treating the molecule with the environmental effects included implicitly. In such reduced description perspectives, information about the molecule will be irreversibly lost to the environment, and the state of the molecule will be accompanied by some degree of uncertainty. Hence, quantum statistical mechanics becomes mandatory to make predictions about the molecule and its properties.? The standard formulation of quantum statistical mechanics is based on (reduced) density operators derived under the approximation of an additively separable composite Hamiltonian, which means that the interaction between the molecule and the environment has been neglected. ?−? ? However, this is a crude approximation, since this includes the terms responsible for the sharing of electrons between the molecule and the environment. Therefore, we will in this work derive a reduced density operator for electronically open molecules, i.e., molecules able to share electrons with the environment, by explicitly averaging out the environmental degrees of freedom of the composite Hamiltonian. This new reduced density operator can be viewed as a generalization of the grand canonical density operator.

In the quantum chemistry community, the standard practice for modeling energetically open molecules, i.e., molecules interacting with the environment through electrostatic interactions, is to incorporate the effect of the environmental interaction into the model of the molecule by means of effective Hamiltonians. ?−? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? It is not yet standard practice to consider electronically open molecules. In response to that, we have recently presented wave function-based ?−? ? and density-operator-based? formalisms targeting electronically open molecules by combining standard electronic structure theory? with open quantum system theory. ?−? ? ? Meanwhile, the term “reduced” appears in different contexts in electronic structure theory and open quantum system theory. Specifically, reduced density matrices (RDMs) are in electronic structure theory routinely used to rewrite expensive many-electron integrals in terms of simpler few-electron integrals to reduce the computational cost of a calculation on a single closed and isolated molecule.? The RDMs may be obtained by integrating out a number of degrees of freedom of the $[eqn]$ -electron density matrix under the condition of $[eqn]$ -representability (i.e., the RDMs must be representable of the full $[eqn]$ -electron problem). ?,? In open quantum system theory, on the other hand, the aim is to describe only a limited number of degrees of freedom (specifically those of the molecule) among a multitude of degrees of freedom (from the environment), and such descriptors are called “reduced”. We target such a reduced descriptor for electronically open molecules in the present work.

The density operator provides the natural framework for quantum statistical mechanics. Specifically, the uncertainty in a molecular system can be implanted in its state description as a statistical mixture (i.e., probabilistically weighted combination) of maximum-possible-information states (typically called pure states), representing the members of the relevant ensemble. ?,? It is customary to distinguish between equilibrium and non-equilibrium quantum statistical mechanics. In equilibrium theory, the weights of the ensemble states are obtained through a maximization of the average entropy subject to constraints related to energy and/or particle conservation. ?−? ? The canonical ensemble permits interfacial energy exchange provided that the average energy of the system stays constant. Analogously, the grand canonical ensemble permits interfacial energy and particle fluctuations given that the average energy and average particle number in the system remain constant. ?−? ?

The vantage point for our derivation of the reduced density operator will be the canonical density operator using the composite Hamiltonian. Specifically, we include the interaction term in the Hamiltonian, which is typically neglected. In order to average out the environmental degrees of freedom, we first have to define the partial trace operation in the composite fermionic Fock space. However, this is a nontrivial task due to the correlation between the molecule’s and the environment’s electronic degrees of freedom, which makes the composite density operator nonseparable and the partial trace operation ambiguous. This is known as the fermionic partial trace ambiguity. ?,? We resolve this ambiguity by working with composite states in a second quantization framework based on a common orthonormal set of orbitals, which allows us to propose an unambiguous definition of the partial trace operation in the composite fermionic Fock space.

In a local spin orbital basis, where the full spin orbital space has been spatially localized onto either the molecule or the environment, the full composite Hamiltonian naturally partitions into a molecule Hamiltonian, an environment Hamiltonian, and an interaction Hamiltonian, which is the structure normally encountered in open quantum system theory. ?−? ? ? Neither the molecule nor the environment Hamiltonian commutes with the interaction Hamiltonian. Therefore, the expansion of the canonical composite density operator will follow from the Zassenhaus expansion, that being an infinite product of exponentiated and increasingly nested commutators. ?,? We approximate commutativity between the subsystem Hamiltonians and the interaction Hamiltonian, such that all commutators in the Zassenhaus expansion vanish. Thereby, the canonical density operator gets effectively partitioned into the components of the open quantum system Hamiltonian. This approximation may be gradually lifted by including more terms in the Zassenhaus expansion, and our methodology thus provides a hierarchy of systematically improvable approximations to the reduced density operator.

By comparing our reduced density operator to the standard grand canonical density operator, we are prompted to define the generalized chemical potential. The naming of this quantity is motivated by its analogous role and its reduction to the chemical potential under the approximations that render our reduced density operator identical to the grand canonical density operator. These approximations are (i) restriction of excitations to occur within the same orbitals and (ii) assumption of equal interaction with the environment for all molecule spin orbitals. The restriction of excitations amounts to reducing all one-electron excitation operators to the corresponding molecule and environment number operators.

The approximation of equal interaction with the environment for all molecule spin orbitals entails the generalized chemical potential to lose its dependence on specific orbitals in the molecule, and thereby, we recover the standard chemical potential. This approximation corresponds to the wide-band approximation (i.e., constant density of states), which is noted to be neither intended nor suitable for molecules. ?−? ? ? ? ? The chemical potential is in standard derivations interpreted as the Lagrange multiplier enforcing the constant average electron number in the molecular system, ?−? ? thus offering no clues on its underlying approximations or physical soundness. Our work alleviates this problem and provides the missing insights into the chemical potential.

The paper is organized as follows. The theoretical preliminaries for our work are presented in Section. Next, we review how to construct states for the composite system, and use this basis to resolve the fermionic partial trace ambiguity by proposing an unambiguous definition of the partial trace operation in the composite fermionic Fock space in Section. The reduced density operator is then derived in Section. Moreover, we compare it with the grand canonical density operator and discuss the physical significance of the approximations invoked in our work. Lastly, we summarize our findings in Section.

Theoretical Preliminaries

2

In this section, we present the theoretical preliminaries for our work. Specifically, we provide a brief overview of quantum statistical mechanics from an electronic structure perspective to introduce the density operator formalism. We then briefly discuss the canonical density operator, since it will be our vantage point for the derivation of the reduced density operator. Lastly, we write up the molecular Hamiltonian in a spatially local spin orbital basis and review how it may be partitioned for two interacting subsystems.

Quantum Statistical Mechanics

2.1

We consider an ensemble of identical molecular systems, in which each of the members is described by the time-independent Schrödinger equation,

$$[eqn]$$

where $[eqn]$ is the standard molecular Hamiltonian (specified in eq), $[eqn]$ is the $[eqn]$ th eigenstate of the molecular system, and $[eqn]$ is the corresponding energy. The molecular system can be described by the density operator defined as ?−? ?

$$[eqn]$$

where $[eqn]$ are the ensemble weights, i.e., the probabilities of the molecular system being in specific eigenstates. The density operator is a convex combination of the eigenstates of the molecular system, meaning that $[eqn]$ and $[eqn]$ . Moreover, we note that eq implicitly assumes that the molecular system is at statistical equilibrium, which means that the corresponding density operator is a stationary solution to the Liouville–von Neumann equation (in atomic units),

$$[eqn]$$

In other words, $[eqn]$ , which implies that $[eqn]$ and $[eqn]$ have simultaneous eigenfunctions. The expectation value of the general operator $[eqn]$ of the molecular system can therefore be computed by?

$$[eqn]$$

The ensemble weights from eq are chosen to maximize the entropy. This ensures that only the available information about the molecular system is used to assign the ensemble weights. The entropy can be represented by the entropy operator,

$$[eqn]$$

and the average entropy can thus be computed by?

$$[eqn]$$

Canonical Density Operator

2.2

A molecular system is rarely alone in space. Therefore, we consider each (large) molecular system as energetically equilibrated (i.e., energetically open to its surroundings and at statistical equilibrium). We can make a molecular system energetically open by allowing for small energy fluctuations between the molecular system and its surroundings, and we can ensure that it remains at statistical equilibrium by enforcing the condition of a constant average energy of the molecular system. This entails a canonical distribution of the molecular systems, which means that the members in the ensemble may have different energies, all within the energy window determined by the size of the allowed energy fluctuations. Maximization of the average entropy over this ensemble results in the ensemble weights $[eqn]$ and thus the canonical density operator, ?−? ?

$$[eqn]$$

where the parameter $[eqn]$ is defined as the derivative of the average entropy of the surroundings at total energy with respect to the energy of the surroundings, ?,? and $[eqn]$ is the canonical partition function, which ensures normalization of the canonical density operator. In the macroscopic (thermodynamic) limit, $[eqn]$ may be interpreted as the inverse temperature (i.e., $[eqn]$ with $[eqn]$ being the Boltzmann constant) due to the definition of temperature.?

The ensemble weights and therefore the form of the canonical density operator in eq rely on the underlying approximation of an additively separable Hamiltonian (i.e., the interaction between the molecular system and its surroundings is neglected). ?−? ? However, we are interested in electronically open molecules, for which we cannot justify such an approximation. As previously mentioned, our vantage point will be a canonically distributed molecular system (and thus the canonical density operator), which we now partition into two interacting subsystems, namely, a small region of interest, henceforth referred to as the molecule, and a large environment surrounding the molecule. The situation is depicted in Figure. The approximation of an additively separable Hamiltonian for the energetically equilibrated molecular system and its surroundings can be justified from our viewpoint (the molecule) because of the nonlocality of the neglected interactions (i.e., the neglected interactions between the molecular system and its surroundings occur far away from the molecule). Hence, our vantage point is physically sound and valid for the purpose of this work.

Illustration of our partitioning of a canonical molecular system energetically equilibrated with its surroundings into two interacting subsystems, namely, the molecule and the environment, which can exchange electrons.

Partitioning the Molecular Hamiltonian

2.3

For a molecular system at a parametrically fixed nuclear geometry (in the Born–Oppenheimer approximation), the Hamiltonian in eq is the standard molecular Hamiltonian. In the spin orbital basis, it takes the following form (excluding nuclear repulsion) in second quantization,?

$$[eqn]$$

The one- and two-electron integrals are defined as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the combined spin and spatial $[eqn]$ coordinate of the $[eqn]$ th electron, $[eqn]$ represents the position vector of the $[eqn]$ th nucleus with nuclear charge $[eqn]$ , and $[eqn]$ are the spin orbitals. The operators $[eqn]$ and $[eqn]$ are the one- and two-electron excitation operators, and $[eqn]$ and $[eqn]$ are the fermionic creation and annihilation operators, respectively. The creation and annihilation operators obey the anticommutation rules for fermionic operators,

$$[eqn]$$ $$[eqn]$$

The spin orbital basis for the composite system $[eqn]$ may now be spatially localized such that the spin orbitals are localized to either the molecule or the environment. Hence, in the local representation, we may write

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the sets of spin orbitals localized to the molecule and the environment, respectively. We note that the molecule and environment spin orbitals are mutually orthogonal, since $[eqn]$ constitutes an orthonormalized set of orbitals. In the absence of covalent molecule-environment bonds, the subsystems have no shared orbitals, and the partitioning in eq is therefore possible. The orthogonality of the orbitals requires that any orbital centered on one subsystem has orthogonalization tails on the other subsystem. Such orthogonalization tails will become detrimental if one chooses basis sets with near-linear dependence, as it results in a nonlocal inverse atomic orbital overlap matrix, $[eqn]$ .? This nonlocality cannot be circumvented as it destroys locality through the presence of $[eqn]$ in a non-orthogonal formulation, or smears the locality of the orbitals if an orthogonal representation is chosen. Consequences of this have been observed, e.g., in transport calculations,? and this has been discussed and resolved, albeit from a different vantage point, by Reuter and Harrison.? Keeping these issues of locality in mind, the approach of orbital space partitioning considered here fully counts an orbital toward the subsystem at which it is centered. Hence, all orbital tails are considered mathematical artifacts? while in reality, small components may be due to constraints of the electronic density matrix (physical requirements). In other words, the size of the tails provides a lower limit on how accurately one may quantify charge-transfer between subsystems in such a framework. The above discussion relies on knowledge acquired ?,? using modern localization functions, which can generate compact occupied and virtual orbitals. ?,?

In the local spin orbital basis in eq, each term in the Hamiltonian for the composite system (eq) can be categorized into belonging to the Hamiltonian for the molecule $[eqn]$ , the Hamiltonian for the environment $[eqn]$ , or the interaction term $[eqn]$ . The full composite Hamiltonian reads

$$[eqn]$$

The terms in $[eqn]$ and $[eqn]$ include only summation over spin orbitals localized to the molecule and environment, respectively,

$$[eqn]$$ $$[eqn]$$

The terms in $[eqn]$ are those that contain mixed summations, and it can be written compactly as

$$[eqn]$$

The interaction Hamiltonian describes the correlation between electrons in the molecule and the environment, and the model choice of both the molecule and environment determines what interactions can be modeled. We have recently presented a categorization of the interaction terms into particle-number conserving and particle-number nonconserving contributions.? For brevity, these terms are referred to as particle-conserving and particle-breaking, respectively. The particle-conserving interactions are responsible for the electrostatic interaction between the subsystems and do not cause fluctuations of electrons between the molecule and the environment. These interactions are routinely considered in the quantum chemistry community through embedding and multiscale schemes. ?−? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? A recent example include the open quantum system theory for polarizable continuum models developed by Guido et al. that include the polarization and dispersion interactions between a solute and solvent.? The particle-breaking components of $[eqn]$ can induce fractional charging of the molecule by connecting Fock space states with different particle numbers. The focus of this work is the construction of an electronically open reduced density operator. Therefore, we omit the particle-conserving interactions in the following, but note that they may be included by means of standard practices.

Recently, we have derived an effective one-electron interaction operator that approximates the particle-breaking interactions.? This interaction operator plays a fundamental role in the present work, and we therefore re-derive and motivate its form in the following. In the local orbital basis (eq), the particle-breaking components of the interaction Hamiltonian may be written as

$$[eqn]$$

The integral $[eqn]$ and its Hermitian conjugate amounts to the overlap between the charge distributions $[eqn]$ and $[eqn]$ , and since the orbitals are separated in space, the charge distributions are vanishingly small. The corresponding overlaps will therefore be negligible, and the particle-breaking interaction Hamiltonian may be approximated as

$$[eqn]$$

where we have anticommuted $[eqn]$ and $[eqn]$ to the right and left of the parentheses, respectively. Lastly, we may wrap up all the integrals by defining the electronic coupling elements,

$$[eqn]$$

with $[eqn]$ being some wave function. The coupling elements represent the strength of the orbital-specific molecule-environment interactions. Thereby, we obtain the following effective one-electron particle-breaking interaction Hamiltonian,?

$$[eqn]$$

In the following, we will refer to the effective one-electron particle-breaking interaction Hamiltonian in eq simply as the interaction Hamiltonian. We note that this bilinear form of the interaction Hamiltonian is typically used in descriptions of bipartite open quantum systems. ?,?−? ? ? ? ? ? ? ? ? ? ? ? ? However, our derivation in ref ? provides a prescription for how to compute the strengths $[eqn]$ in terms of standard one- and two-electron integrals. On the downside, the computation of the coupling elements requires a full quantum mechanical description of the composite system, and we are not interested in performing such calculations. Therefore, we have in our previous works ?−? ? ? estimated the coupling elements based on the expected magnitudes of the contained interactions. In the future, we will be able to compute and benchmark these coupling elements for small composite systems using the charge-localized electronic wave functions recently developed by Folkestad and Høyvik.?

Furthermore, we note that there is a direct connection between the presented formalism and electron transfer theory. In the seminal Marcus electron transfer theory,? the electron transfer rates are related to the electronic coupling elements between an acceptor and a donor (a redox pair), which are both treated explicitly. ?−? ? Central to the calculation of these electronic coupling elements is the Hamiltonian matrix element, which couples two diabatic (or rather, quasi-diabatic) states. ?−? ? ? In the so-called direct methods, these coupling elements are evaluated using states obtained by diabatization of the corresponding adiabatic states ?−? ? ? ? ? ? or computed directly, e.g., from constrained density functional theory. ?−? ? ? ? The concepts used to develop the theory for electronically open molecules may similarly be used to calculate such coupling elements, since the formalism is based on the construction of quasi-diabatic states for the molecule. Hence, for the application in electron transfer theory, one may construct quasi-diabatic states for the donor and the acceptor, rather than considering a molecule and the environment. Thereby, the particle-breaking terms of $[eqn]$ (in this work approximated by $[eqn]$ ) become responsible for the electronic coupling elements used in electron transfer theory. This is discussed in more detail by Folkestad and Høyvik in ref ?.

Composite Fermionic Fock Space

3

A prerequisite for our work is the ability to average over the environmental degrees of freedom in the composite fermionic Fock space. However, the requirement of anticommutation between electrons in the molecule and the environment makes this nontrivial. The composite fermionic Fock space is an antisymmetric tensor product (normally called a wedge product) of the molecule and environment Fock spaces. The structure of the wedge product ensures the anticommutation between electrons in the molecule and the environment, but by partially tracing out the environmental degrees of freedom of the composite density operator, the information necessary to enforce the global anticommutation rules (eqs and ?) is lost. This problem is known as the fermionic partial trace ambiguity. ?,?

Szalay et al. have recently presented an approach to integrate the anticommutativity of fermionic modes between disjoint sets into the individual matrix elements of the composite density operator, such that the antisymmetric structure of the composite fermionic space is preserved when partially tracing out the noninteresting subsystems.? Specifically, they apply a Jordan–Wigner transformation? to effectively strip the electrons of their fermionic behavior. The result is a bosonic-like structure containing the information about the fermionic occupations and matrix element-specific sign factors ensuring the fulfillment of the global anticommutation rules. The fermionic nature of the resulting reduced density operator may then be restored through the inverse transformation.?

In the present work, we provide an alternative strategy based on the definition of composite states in a second quantization framework built from a common orthonormal set of orbitals. This setup allows us to unambiguously define the partial trace operation in the composite fermionic Fock space, such that the information required to enforce the global anticommutation rules is absorbed into the matrix elements of the reduced density operator. Thereby, we resolve the fermionic partial trace ambiguity.

This section is organized as follows. First, we review how to construct a composite basis for the interacting composite system from the eigenstates of the noninteracting Hamiltonians. This composite basis is then used to define the fermionic partial trace operation.

States for the Noninteracting Molecule and

Environment

3.1

The isolated molecule and environment Hamiltonians describe the noninteracting molecule and environment, respectively, and their eigenstates can be written as

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ $[eqn]$ and $[eqn]$ $[eqn]$ being an eigenstate and corresponding energy eigenvalue of the molecule (environment), e.g., the full configuration interaction eigenpair.

We may further define the molecule and environment number operators as

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ and $[eqn]$ are the corresponding one-electron number operators, and we note that $[eqn]$ and $[eqn]$ are orthonormal subsets of the orthonormal composite basis $[eqn]$ (eq). The total number operator is the sum of the molecule and environment number operators,

$$[eqn]$$

It can be shown that $[eqn]$ and $[eqn]$ commute with their respective number operators,

$$[eqn]$$

Hence, all eigenstates $[eqn]$ and $[eqn]$ are eigenstates of $[eqn]$ and N̂ E, respectively, and the states may be additionally labeled according to the number of electrons in the state,

$$[eqn]$$ $$[eqn]$$

In the following, the electron numbers are omitted in the state labels to simplify the notation, but we include eqs and ? in order to (i) highlight that each state in the sets $[eqn]$ and $[eqn]$ corresponds to a particle number eigenstate and (ii) emphasize that we in this work consider all Fock states, not just those that belong to a particular number of electrons.

Many-Electron Basis for the Composite System

3.2

The molecule and environment Hamiltonians commute $[eqn]$ . Therefore, simultaneous eigenstates may be constructed for $[eqn]$ from the sets $[eqn]$ and $[eqn]$ . For that purpose, we define state operator strings for the molecule $[eqn]$ and the environment $[eqn]$ ,

$$[eqn]$$

where we use the notation $[eqn]$ and $[eqn]$ to denote the spin orbitals occupied in the particular molecule and environment state, respectively. The composite states may be defined by means of the state operator strings as^47^

$$[eqn]$$

with the Hermitian adjoint of the state defined by

$$[eqn]$$

The order in which the molecule and environment state operator strings act during the preparation of the composite state is one choice, and we note that all other arrangements would have been equally valid. This leads to a sign ambiguity (called the fermionic ambiguity) in composite states where only one of the subsystems contains an odd number of electrons because of the anticommutation rules in eq. ?,? However, the composite density operator is unaffected by the fermionic ambiguity, since its construction eliminates all potential sign ambiguities. Moreover, we note that the orbitals from which the state operator strings $[eqn]$ and $[eqn]$ are constructed belong to the overall orthonormal set in eq. Therefore, the set $[eqn]$ defines an orthonormal many-electron basis. The composite states $[eqn]$ are eigenstates of $[eqn]$ ,

$$[eqn]$$

but they are not eigenstates of the full or approximate composite Hamiltonian ( $[eqn]$ or $[eqn]$ ). The composite states will, in general, not be eigenstates of the total spin operator $[eqn]$ and will therefore lack spin symmetry. However, if needed, we may construct configuration state functions (CSFs) (i.e., spin-adapted combinations) to enforce particular spin symmetries. The wave function of the composite system may then be expanded in either the composite states $[eqn]$ or the basis of the CSFs. In the present work, the composite wave function will not be constructed explicitly, as we will only use its existence in the derivation of the reduced density operator.

Fermionic Partial Trace Operation

3.3

The density operator for the composite system may be written in the basis of the composite states defined in eqs and ?,?

$$[eqn]$$

The matrix elements $[eqn]$ are nonseparable due to the correlation between the electronic degrees of freedom in the molecule and the fermionic environment. In other words, the composite density operator cannot be written as a normal tensor product of fermionic density operators representing the molecule and the environment. ?,?,?

Nonetheless, the troublesome matrix elements $[eqn]$ in eq are typically approximated as being separable such that the reduced density operator of the molecule can be obtained from the composite density operator by averaging over (i.e., tracing out) the environmental degrees of freedom. ?,?−? ? ? ? ? ? ? The accuracy of this approximation can be checked by performing a Schmidt decomposition of the composite state, which may be considered approximately separable if only one Schmidt coefficient is dominant.? The advantage of our composite basis (Section) is that we can unambiguously define the partial trace operation. Specifically, we propose the following definition of the partial trace operation in the composite fermionic Fock space,

$$[eqn]$$

We note that this definition can similarly be applied without loss of generality to operators defined on composite bosonic-bosonic and fermionic-bosonic tensor product spaces.

Our partial trace operation may now be used to trace out the environment states of the composite density operator in eq. The result is the reduced density operator for the molecule,

$$[eqn]$$

with

$$[eqn]$$

The problem of decomposing the nonseparable matrix elements $[eqn]$ , such that the antisymmetric structure of the composite space is preserved, has just been turned into how one can construct and compute the reduced matrix elements in eq. We can build and compute (or rather, approximate) such matrix elements using our partitioning of the Hamiltonian (Section) and definition of composite states (Section) combined with the framework of standard electronic structure theory. As previously mentioned, the fermionic partial trace ambiguity arises when the information required to enforce the anticommutation rules between electrons in the molecule and the environment is lost. This information is contained in the reduced density matrix elements. Consequently, we argue that our definition of the partial trace operation (eq) together with our definition of composite states (Section) resolves the fermionic partial trace ambiguity.

Brief Note on Covalent Interactions

3.3.1

The existence of the spatially localized spin orbital basis in eq is guaranteed as long as the molecule and environment are not connected by a covalent bond (i.e., share orbitals). Therefore, our previous work ?−? ? ? has focused on describing noncovalent equilibrium interactions, e.g., molecules interacting with solvent molecules or physisorbed on molecular surfaces, but we plan to extend our formalism to covalent bonding situations in future work. However, we note that the presence of covalent bonds will introduce another ambiguity in the reduced description of the molecule related to the partial trace operation. In the following, we comment on this ambiguity in order to (i) avoid confusion with the fermionic partial trace ambiguity (i.e., loss of anticommutativity between electrons in the molecule and the environment) and (ii) emphasize that we in this work completely avoid this ambiguity by requiring the interaction between the molecule and the environment to be noncovalent.

The ambiguity arises when distinct composite states give rise to the same reduced description. ?,? For example, a bonding interaction between the molecule and an environment may be described as a linear combination of Slater determinants built from orbitals which are strictly localized to the molecule or environment,

$$[eqn]$$

The symmetric and antisymmetric combination represents ground and excited states. The composite wave function is a pure state, and the composite density operators can be written as

$$[eqn]$$

where we have used the notation $[eqn]$ and $[eqn]$ . However, by partially tracing out the environmental degrees of freedom of each composite density operator, we obtain the same reduced description,

$$[eqn]$$

and the information detailing the character of the state is lost. This creates an ambiguity in the reduced description, and we note that it is not unique for fermionic-fermionic bipartite systems, but also applies to bosonic-bosonic and fermionic-bosonic composite systems.

Reduced Density Operator

4

In this section, we derive the reduced density operator for electronically open molecules by explicitly averaging out the environmental degrees of freedom of the canonical density operator (eq) equipped with the approximate composite Hamiltonian $[eqn]$ . The derivation is provided in Section, and we construct the corresponding reduced partition function in Section. We then identify the approximations that render our reduced density operator identical to the grand canonical density operator in Section. By comparing our reduced density operator to the grand canonical density operator, we are prompted to define a quantity that we coin the generalized chemical potential, and its physical significance is elaborated upon in Section. Our derivation of the reduced density operator relies on an approximate commutativity between the subsystem Hamiltonians and the interaction Hamiltonian, and we discuss this assumption in Section.

Derivation

4.1

We now derive the reduced density operator for an electronically open molecule. Our vantage point will be the canonical composite density operator, i.e., the canonical density operator (eq) with the approximate $[eqn]$ composite Hamiltonian (eq),

$$[eqn]$$

The expansion of the canonical composite density operator follows from the Zassenhaus expansion. ?,? Meanwhile, to ensure that the resulting reduced density operator obtains a Hermitian form, we write the exponential as

$$[eqn]$$

The first exponential in the second equality of eq may similarly be expanded,

$$[eqn]$$

$[eqn]$ , $[eqn]$ , and $[eqn]$ are all Hermitian operators, and we may consider the Hermitian adjoint,

$$[eqn]$$

where we have used that

$$[eqn]$$

to write the last equality. We note that

$$[eqn]$$

Hence, we may substitute the first exponential in eq with eq and obtain

$$[eqn]$$

The subsystem Hamiltonians are now approximated to commute with the interaction Hamiltonian, an approximation that will be discussed in more detail in Section,

$$[eqn]$$

Thereby, all commutators vanish, and the composite density operator is effectively split into the components of the composite Hamiltonian,

$$[eqn]$$

where we note that $[eqn]$ since $[eqn]$ . $[eqn]$ is the approximate canonical partition function whose role is to ensure normalization of our approximate canonical composite density operator,

$$[eqn]$$

To obtain the reduced density operator, we use the definition of the partial trace operation introduced in Section. Partially tracing out the environmental degrees of freedom of the approximate canonical composite density operator in eq yields

$$[eqn]$$

where we have used that $[eqn]$ and $[eqn]$ are eigenstates of $[eqn]$ and $[eqn]$ , respectively. The exponential of the interaction Hamiltonian is expanded in its Taylor series. Hence, the transition moment reads

$$[eqn]$$

In the following, we evaluate $[eqn]$ and $[eqn]$ . The first transition moment becomes

$$[eqn]$$

where we have realized that the composite states contain identical environment states and that the interaction Hamiltonian will either remove or add an electron to the environment. The action of the interaction Hamiltonian will thus result in a vanishing overlap between environment states of different particle-number symmetry. In other words, there is no first-order contribution to the reduced density operator, which is consistent with typical open quantum system treatments. ?−? ? ? The second transition moment reads

$$[eqn]$$

Next, the resulting transition moments in eq are evaluated. The first one becomes

$$[eqn]$$

In the first equality, we have used the definition of the composite states (eq) and moved $[eqn]$ to the left of $[eqn]$ . In the second equality, we have similarly reordered the state operator strings by molecule and environment degrees of freedom. The reordering of the state operator strings results in the overall sign $[eqn]$ . Meanwhile, we note that state $[eqn]$ and $[eqn]$ must contain the same number of electrons (i.e., $[eqn]$ ); otherwise, the overlap will be zero. The combined sign factor thus reduces to $[eqn]$ . In addition, we have inserted the resolution of identity in the composite Fock space (i.e., $[eqn]$ ) in the third equality. By using the definition of the composite states and noting that $[eqn]$ and $[eqn]$ must both be the identity for the overlap to be nonzero, we obtain the fourth equality. Lastly, we let the state operator strings work on the vacuum states, thereby separating the expectation value into a molecule and environment part. In a similar way, the second transition moment in eq evaluates to

$$[eqn]$$

To simplify eq by means of the results of eqs and ?, we first define the canonical-like density operator of the environment,

$$[eqn]$$

with the environment partition function,

$$[eqn]$$

Thereby, we may write eq as

$$[eqn]$$

We now identify the environment average,

$$[eqn]$$

Thereby, eq becomes

$$[eqn]$$

Next, we anticommute $[eqn]$ and $[eqn]$ , such that we can simplify the reduced density operator expression by defining a generalized chemical potential,

$$[eqn]$$

and the molecule orbital energy shift,

$$[eqn]$$

The naming of $[eqn]$ is elaborated upon in Section. In terms of the generalized chemical potential (eq) and the molecule orbital energy shift (eq), the reduced density operator reads

$$[eqn]$$

The terms between $[eqn]$ and $[eqn]$ resemble the first two terms of the Taylor series of an exponential, and we may approximately recollect the terms in an exponential accordingly. The second-order term of this Taylor series expansion will arise from $[eqn]$ , since all odd powers of $[eqn]$ result in overlaps between environment states of different particle-number symmetry, which vanish as argued in eq. However, we show in Section that the first nonvanishing correction to our approximation of commutativity between the subsystem Hamiltonians and the interaction Hamiltonian (eq) originates from an effective $[eqn]$ -type of interaction. Therefore, we decide not to consider any higher-order terms of the Taylor series expansion in eq to stay consistent with this approximation.

In addition to collecting the terms in an exponential, we now use the definition of the resolution of the identity (in the molecule Fock space) to bring the reduced density operator into the form,

$$[eqn]$$

where we have also identified the reduced partition function for the molecule, $[eqn]$ . This reduced partition function is considered in more detail in Section. We note that the orbital energy shift (defined in eq) is a constant scalar and commutes with all operators. It may therefore be moved and absorbed into the molecule Hamiltonian,

$$[eqn]$$

The molecule orbital energy shift directly results from the relaxation induced by the coupling of the molecule to the interacting environment. Consequently, the reduced density operator becomes

$$[eqn]$$

where we have defined the generalized chemical potential operator,

$$[eqn]$$

The molecule Hamiltonian is total symmetric (i.e., spin-free) and $[eqn]$ may similarly be constrained to be of singlet symmetry in a spin-restricted basis. In this way, the reduced density operator can be restricted to be of singlet spin symmetry when desirable, and hence only connect states of the same spin symmetry.

Reduced Partition Function

4.2

The reduced partition function $[eqn]$ ensures the normalization of the reduced density operator. The definition of $[eqn]$ is given in eq, and the canonical composite partition function defined in eq can be written as

$$[eqn]$$

where we have exploited that $[eqn]$ and $[eqn]$ are eigenstates of $[eqn]$ and $[eqn]$ , respectively. The expectation value in eq is evaluated following the same procedure as shown from eq to ?, the only difference being $[eqn]$ instead of $[eqn]$ , which directly entails the overall sign $[eqn]$ . The reduced partition function therefore becomes

$$[eqn]$$

where the third equality follows from absorbing the constant molecule orbital energy shifts into the molecule Hamiltonian, as done in eq.

Grand Canonical Density Operator

4.3

We now want to connect our reduced density operator (eq) to the grand canonical density operator. In other words, we want to identify the approximations that will render our reduced density operator identical to the grand canonical density operator. This will provide insights into the underlying approximations and physical content of the grand canonical density operator and the chemical potential.

By restricting the excitation operators in eq to occur within the same orbitals, the one-electron excitation operators for the molecule and environment reduce to the corresponding molecule and environment one-electron number operators as defined in eqs and ?, respectively. Specifically, the restriction of the environment one-electron excitation operators $[eqn]$ entails the generalized chemical potential to read

$$[eqn]$$

The restriction on the molecule one-electron excitation operators $[eqn]$ results in

$$[eqn]$$

$[eqn]$ are the diagonal elements of eq and can be written as

$$[eqn]$$

where we have introduced the one-electron chemical potential,

$$[eqn]$$

The one-electron chemical potential contains information on the flow of electrons in non-equilibrium situations and will be discussed in Section. Similarly, under this restriction on the excitation operators, the molecule orbital energy shift $[eqn]$ used in the redefinition of the molecule Hamiltonian becomes

$$[eqn]$$

Finally, by approximating all molecule spin orbitals to have equal interaction with the environment (i.e., $[eqn]$ ), the generalized chemical potential is the same for all molecule spin orbitals, and we may write

$$[eqn]$$

This approximation amounts to the wide-band approximation known from non-equilibrium transport theory. ?−? ? ? Altogether, the approximated reduced density operator becomes

$$[eqn]$$

where we have also observed that the reduced partition function reduces to the standard grand canonical partition function under these approximations (i.e., $[eqn]$ ). In addition, we have used that the molecule Hamiltonian and molecule number operator commute (eq) to gather the terms in a single exponential in the last equality. In contrast, $[eqn]$ , which means that the reduced density operator in eq cannot be collected in one exponential. The result of eq is the grand canonical density operator.

We note that expectation values of the environment one-electron excitation operators naturally reduce to the environment number operators in the case of a single-determinantal environment wave function. However, this is not the case for the molecule one-electron excitation operators nor if multideterminantal wave functions, such as coupled cluster and configuration interaction wave functions, are used for describing the environment.

Generalized Chemical Potential

4.4

An important feature of our reduced density operator (eq) relative to the grand canonical density operator (eq) is the independence of the generalized chemical potential with respect to the level of theory used to describe the environment. In other words, our reduced density operator is fully compatible with all possible level-of-theory treatments of the environment, ranging from full configuration interaction theory to Fermi–Dirac statistics. ?,? This is in stark contrast to standard protocols, where the chemical potential determines the environment orbital occupancy, ?−? ? and not the other way around as in our framework.

The diagonal values of the generalized chemical potential (eq) are interpreted as the energy released or absorbed by the molecule upon the removal or addition, respectively, of a (fractional) electron. We denote this process as direct electron transfer. This aligns with the standard interpretation of the chemical potential, apart from the possibility of fractional rather than strictly integer electron transfer in our framework. ?−? ? The off-diagonal elements are interpreted as the energy contribution following the coupling to the environment on the excitations within the molecule. This process is termed indirect electron transfer.

Insight into the direction of electron flow in non-equilibrium situations can be obtained from our definition of the generalized chemical potential. Specifically, a positive chemical potential (i.e., energy is released) represents the situation where an electron can be transferred from the molecule to the environment, while a negative chemical potential (i.e., energy is absorbed) corresponds to an electron being transferred from the environment to the molecule. In other words, the direction of the electron flow in non-equilibrium situations is dictated by the sign of the chemical potential, as is standard following the principle of electronegativity equalization. ?,? The sign of our generalized chemical potential directly follows from the value of $[eqn]$ . The average orbital occupancy of the environment can take any value in the range $[eqn]$ , which implies that $[eqn]$ . Consequently, if the environment orbital on average is unfilled or less than half-filled, the generalized chemical potential will be positive, and electron transfer from the molecule to the environment is favored [Figure(left)]. On the other hand, if the environment orbital on average is more than half-filled or filled, the generalized chemical potential will be negative, and electron transfer from the environment to the molecule is favored [Figure(right)]. A special situation arises when $[eqn]$ , entailing $[eqn]$ , which represents an equilibrium situation between electron addition and removal [Figure(center)]. This aligns with our expectations, and the derivation reveals that the direction of electron flow is built directly into our definition of the generalized chemical potential.

Direction of electron flow in non-equilibrium situations can be predicted by the average environment orbital occupation. The different situations are favored electron transfer to the environment (left), equilibrium between electron addition and removal (center), and favored electron transfer to the molecule (right). The molecule’s ability to donate or accommodate an electron is taken into account by the molecule number operator. Therefore, the direction of electron flow is independent of the specific molecule spin orbitals.

Lastly, we note that the equation for the diagonal elements of our generalized chemical potential (i.e., direct electron transfer, eq) resembles the Fermi golden rule? type of equation typically encountered in electron transfer theory. ?−? ? Specifically, it depends on the squared coupling elements and an internal occupation number condition to ensure a physically meaningful transfer process, i.e., whether the environment can donate or accommodate the transferred electron. The molecule’s ability to donate or accommodate an electron is taken into account by the molecule number operator, and we note that overall particle conservation is built into the equations. This may be demonstrated by writing the diagonal element of the first nonvanishing coupling element responsible for the electron transfer (eq) as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the electron numbers in orbitals $[eqn]$ and $[eqn]$ in states $[eqn]$ and $[eqn]$ , respectively. Hence, eq shows that direct electron transfer can only occur if the molecule and environment spin orbital pair constitute one occupied and one virtual spin orbital. This is the particle-number equivalent of the resonance (i.e., energy conservation) condition in Fermi’s golden rule.

Noncommutativity Effects

4.5

The subsystem Hamiltonians were approximated to commute with the interaction Hamiltonian in eq. This may be a strong assumption, and careful numerical studies are required to investigate the validity of this approximation. However, we note that this approximation is implicitly adopted in standard derivations of the grand canonical density operator, and that our work opens up the opportunity to study this assumption. Specifically, our methodology provides a hierarchical approach to improving this approximation (and hereby the reduced density operator) through the inclusion of higher-order (with respect to $[eqn]$ ) terms in the expansion. In the limit of the full Zassenhaus expansion, the expansion is exact but intractable to work with. Therefore, we need approximations, and the assumption of commutativity invoked in Section amounts to the lowest rung.

We may perform a preliminary investigation by identifying the first correction to the noncommutativity between the subsystem Hamiltonians and the interaction Hamiltonian. For that purpose, we retain the next higher-order terms in eq, such that the approximate canonical composite density operator now reads

$$[eqn]$$

with $[eqn]$ being the corresponding approximate canonical composite partition function. Thereby, the reduced density operator becomes

$$[eqn]$$

The exponentials within the transition moments are expanded in their respective Taylor series and regrouped,

$$[eqn]$$

A comparison of eq with eqs and ? reveals that the first nonvanishing correction to the reduced density operator is

$$[eqn]$$

The first correction to the noncommutativity between the subsystem Hamiltonians and the interaction Hamiltonian is noted to contain an effective $[eqn]$ -type of interaction (similar to that contained in our reduced density operator in eq) but at one additional power of $[eqn]$ . Numerical studies of the commutators and $[eqn]$ are required to comment on the magnitude of this correction term.

Summary

5

In this work, we derive a reduced density operator for electronically open molecules by explicitly tracing out the environmental degrees of freedom of the composite Hamiltonian. Specifically, we include the particle-number nonconserving (particle-breaking) interactions responsible for electron fluctuations between the molecule and the environment, which are neglected in standard formulations of quantum statistical mechanics. The derivation is based on composite states in a second quantization framework built from a common orthonormal set of orbitals. This allows us to propose an unambiguous definition of the partial trace operation in the composite fermionic Fock space. Thereby, we resolve the fermionic partial trace ambiguity.

The full composite Hamiltonian naturally splits into a molecule Hamiltonian, an environment Hamiltonian, and an interaction Hamiltonian in a spatially localized spin orbital basis. The new reduced density operator is based on the approximation of commutativity between the subsystem Hamiltonians (i.e., molecule and environment Hamiltonians) and the interaction Hamiltonian. This entails the canonical composite density operator to be effectively partitioned into a molecule, an environment, and an interaction component. However, this assumption may be gradually lifted by including more terms in the Zassenhaus expansion, and our methodology thus provides a hierarchy of systematically improvable approximations to the reduced density operator.

The new reduced density operator can be viewed as a generalization of the grand canonical density operator. Specifically, we identify the approximations that render our reduced density operator equal to the grand canonical density operator. These approximations are (i) restriction of excitations to occur within the same orbitals and (ii) assumption of equal interaction with the environment for all molecule spin orbitals (i.e., the wide-band approximation). This provides new insights into the inherent approximations and physical content of the grand canonical density operator.

We are prompted to define the generalized chemical potential. The naming of this quantity is motivated by its analogous role and its reduction to the standard chemical potential under the approximations that connect our reduced density operator to the grand canonical density operator. The diagonal elements of the generalized chemical potential are interpreted as the energy associated with the addition or removal of a (fractional) electron. This aligns with the standard interpretation of the chemical potential, apart from the possibility of fractional rather than strictly integer electron transfer in our framework. The chemical potential is in standard formulations of quantum statistical mechanics interpreted as the Lagrange multiplier enforcing the constant average electron number in the molecular system. In other words, the chemical potential determines the electron occupancy (top-down) in standard approaches. Meanwhile, our (bottom-up) framework enables an explicit consideration of the electron occupancy in the environment at any level of theory, irrespective of the model used to describe the molecule. Consequently, our reduced density operator is fully compatible with all possible level-of-theory treatments of the environment.

We obtain insights into the direction of electron flow in non-equilibrium situations from our definition of the generalized chemical potential. We find that if the environment orbital on average is unfilled or less than half-filled, the generalized chemical potential will be positive, and electron transfer from the molecule to the environment is favored. Similarly, if the environment orbital on average is more than half-filled or filled, the chemical potential will be negative, and electron transfer from the environment to the molecule is favored. If the environment orbital is half-filled on average, the chemical potential is zero, and this represents an equilibrium situation between electron addition and removal. The equation for the diagonal elements of the generalized chemical potential (i.e., direct electron transfer) resembles the Fermi golden rule type of equation, and we show that overall particle conservation is built into the equations. This is the particle-number equivalent of the resonance (i.e., energy conservation) condition in Fermi’s golden rule.

Lastly, we note that our derivation is based on the assumption of a noncovalent equilibrium interaction between the molecule and the environment. It is of great interest to extend this to include covalent bonding and non-equilibrium situations, and we will in future work explore such extensions.