Study of Nonequilibrium Size and Concentration Effects on the Heat and Mass Diffusion of Indistinguishable Particles Using Steepest-Entropy-Ascent Quantum Thermodynamics

Li, Guanchen; von Spakovsky, Michael R.

doi:10.1115/1.4036735

Conventional first-principle approaches for studying nonequilibrium processes depend on the mechanics of individual particles or quantum states and as a result require many details of the mechanical features of the system to arrive at a macroscopic property. In contrast, thermodynamics, which has been successful in the stable equilibrium realm, provides an approach for determining macroscopic properties without the mechanical details. Nonetheless, this phenomenological approach is not generally applicable to a nonequilibrium process except in the near-equilibrium realm and under the local equilibrium and continuum assumptions, both of which limit its ability to describe nonequilibrium phenomena. Furthermore, predicting the thermodynamic features of a nonequilibrium process (of entropy generation) across all scales is difficult. To address these drawbacks, steepest-entropy-ascent quantum thermodynamics (SEAQT) can be used. It provides a first-principle thermodynamic-ensemble based approach applicable to the entire nonequilibrium realm even that far-from-equilibrium and does so with a single kinematics and dynamics, which crosses all temporal and spatial scales. Based on prior developments by the authors, SEAQT is used here to study the heat and mass diffusion of indistinguishable particles. The study focuses on the thermodynamic features of far-from-equilibrium state evolution, which is separated from the specific mechanics of individual particle interactions. Results for nonequilibrium size (volume) and concentration effects on the evolutionary state trajectory are presented for the case of high temperature and low particle concentration, which, however, do not impact the generality of the theory and will in future studies be relaxed.

Introduction

In addition to the phenomenological approach of studying heat and mass diffusion using empirical equations, there also exist numerous methods for modeling such nonequilibrium phenomena from a first-principle standpoint with each restricted to its own applicable set of spatial and temporal scales. At the macroscopic level, continuum nonequilibrium thermodynamics with the local equilibrium assumption is used but cannot generally be applied at atomistic/mesoscopic scales since the small dimensions of the system result in quantum and for that matter classical effects that the continuum assumption cannot address. Moreover, nonequilibrium processes in the far-from-equilibrium realm make the application of the continuum formulation of nonequilibrium thermodynamics, i.e., the so-called Onsager formulation (e.g., see Ref. [1]), questionable due to its underlying assumption of linearity or near-equilibrium behavior. Furthermore, each specific set of spatial and temporal scales usually entails a different kinematic and dynamic description of system state and its motion. Thus, a general approach which could provide a thermodynamic analysis of nonequilibrium evolution, especially that far-from-equilibrium, across different spatial and temporal scale has been lacking, even though the general system properties of energy and entropy are well defined [2] and their evolutions observable.

To arrive at a general study of nonequilibrium state evolution (i.e., of the entropy generation process), a mathematical framework with a single kinematic and dynamic description that crosses all the temporal and spatial scales and accounts for both noncontinuum quantum and classical effects is needed. At the same time, retaining the advantage of the equilibrium thermodynamic approach, which is able to provide system property information phenomenologically, is also desirable as is providing thermodynamic features of system state and state evolution in the sense of a nonequilibrium relaxation pattern (in the sense of GENERIC [3,4]), which captures the complex and detailed dynamic balance of particles or quantum states independent of the exact details of the micromechanical interactions taking place. Such a framework must also avoid the computational burdens inherent to existing methods based on mechanics (e.g., the Boltzmann equation [5] and molecular dynamics [6]) or quantum mechanics (e.g., “open-system” quantum thermodynamics [7]) that require detailed interaction information of particles or quantum states.

A mathematical framework that meets these requirements is that of SEAQT [8–13]. It can be used to model nonequilibrium processes (even those far-from-equilibrium) from the atomistic to the macroscopic level [14–19]. SEAQT bases its framework on properties well defined at all scales for equilibrium as well as nonequilibrium states such as the energy, particle number, and entropy [2]. The nonlinear dynamics of state evolution are characterized by the entropy generation, which results from the principle of steepest entropy ascent (SEA). This principle forms the basis of the equation of motion that tracks the evolution of energy and entropy in state space. Using the concept of hypoequilibrium state (i.e., a nonequilibrium relaxation pattern), the nonequilibrium trajectory of system state evolution for a large range of initial conditions can be fully described [16–18,20]. In this way, the thermodynamic analysis of both reactive and nonreactive nonequilibrium phenomena, including heat and mass diffusion, can be studied within a single framework. This characteristic of SEAQT enables the comparative study of nonequilibrium trajectories for coupled and uncoupled heat and mass diffusion at different spatial scales. In addition, SEAQT is a first-principle, thermodynamic-ensemble based approach, which views the nonequilibrium relaxation process from the point of view of an ensemble entity and as a result avoids the computational burdens inherent to existing mechanical approaches. Its equation of motion is, thus, able to provide trajectories in thermodynamic state space independent of the mechanical details of the relaxation process, i.e., of the details of specific microscopic interactions. Nonequilibrium features of system behavior are nonetheless clearly revealed.

In this paper, how concentration and system size influence the heat and mass diffusion process is illustrated. Note that the SEAQT framework requires no assumptions of equilibrium nor of near-equilibrium. Furthermore, system energy eigenstructure changes and entropy generation play key roles since system dimensions influence the eigenstructure, while the eigenstructure determines the topology of eigenstates that can be accessed by any given thermodynamic state at any instant of time along the nonequilibrium trajectory of state evolution. Size and concentration effects can, thus, be studied from the point of view of thermodynamics. The present study focuses on the effects of concentration and system size on the kinetics of the nonequilibrium state evolution, specifically, the trajectory in thermodynamic state space. The discussion is limited to the condition that the temperature is sufficiently high and the concentration of the particle is not too great so that Bose–Einstein condensation does not occur. These conditions will be relaxed in future studies but do not impact the generality of the theory described here. Results for fermions and bosons are included here. Two kinds of size effects are presented, i.e., those due to concentration and those due to volume.

This paper starts with an introduction of the SEAQT equation of motion in Sec. 2. Next, the important concepts of the kinetics and dynamics of the nonequilibrium evolution are presented in Sec. 2.1. This is followed by an illustration of hypoequilibrium and nonequilibrium intensive properties in order to give a complete description of the nonequilibrium state evolution in Sec. 2.2. Section 3 then provides a discussion of how interactions between systems are captured within the SEAQT framework followed by results from the application of the SEAQT framework to the prediction of the heat and mass diffusion of fermion and boson systems in Sec. 4. Section 5 then ends with some conclusions.

SEAQT Equation of Motion

SEAQT Equation of Motion.

In this section, the system and state description in SEAQT is given, and the equation of motion is presented. According to the discussion by Grmela et al. [3,4,21] and Beretta and coworkers [12,22], a nonequilibrium framework is a combination of both irreversible relaxation and reversible symplectic dynamics in general. The equation of motion takes the generalized form of the Ginzburg–Landau equation [3,22]

\frac{d}{d t} θ (t) = X_{θ (t)}^{H} + Y_{θ (t)}^{H}

(1)

where

θ (t)

represents the state evolution trajectory and

X_{θ (t)}^{H}

and

Y_{θ (t)}^{H}

are the functions of the system state

θ (t)

and represent the reversible symplectic dynamics and irreversible relaxation process, respectively. In the SEAQT framework, the system is defined by the Hamiltonian

\hat{H}

⁠, the system state is represented by the density operator

\hat{ρ}

(or probability distribution in phase space [18]),

X_{θ (t)}^{H}

follows the Schrödinger equation (or Louisville's equation [18]), and

Y_{θ (t)}^{H}

is derived from the SEA principle. To describe the evolutionary process, we use conservation laws to construct the equation of motion, following the approach given by Ref. [23]. The result is

\frac{d \hat{ρ}}{d t} = \frac{1}{i ℏ} [\hat{ρ}, \hat{H}] + \frac{1}{τ} \hat{D}

(2)

where the first term is the Schrödinger term, and the second the dissipation term. If the system is in a pure (zero-entropy) state initially, i.e.,

\hat{ρ} \hat{ρ} = \hat{ρ}

⁠, the equation of motion reverts back to the Schrödinger equation of quantum mechanics. If the system is in a so-called mixed (nonzero-entropy) state initially and

\hat{ρ}

is diagonal in the energy eigenstate basis,

\hat{ρ}

corresponds to a classical probability distribution in energy, and

\hat{H}

commutes with

\hat{ρ}

so that the reversible term goes to zero even though

\hat{ρ}

may not be a Maxwellian distribution among the energy eigenlevels. The state evolution of such a mixed-state operator cannot be captured by the reversible Schrödinger term and is instead given by the second term to the right of the equals, the dissipation term, which captures the probability redistribution toward the Maxwellian distribution. The dissipation term is constructed using a set of operators called the “generators of the motion.” Each generator corresponds to one of the conservation laws to which the system is subjected. For example, an isolated system is subject to two conservation laws, probability normalization and energy conservation, so that the generators of the motion are

{\hat{I}, \hat{H}}

⁠, which are the identity operator (⁠

T r {\hat{I}, \hat{ρ}} = 1

⁠) and the Hamiltonian (⁠

T r {\hat{H}, \hat{ρ}} = energy

⁠). For the case when

\hat{ρ}

is diagonal in the energy eigenstate basis, i.e., only dissipation term is active, the equation of motion takes the form [10]

\frac{d p_{k}}{d t} = \frac{1}{τ} \frac{| \begin{matrix} - p_{k} \ln p_{k} & p_{k} & ε_{k} p_{k} \\ 〈 s 〉 & 1 & 〈 e 〉 \\ 〈 e s 〉 & 〈 e 〉 & 〈 e^{2} 〉 \end{matrix} |}{| \begin{matrix} 1 & 〈 e 〉 \\ 〈 e 〉 & 〈 e^{2} 〉 \end{matrix} |} = \frac{1}{τ} D_{k} (p)

(3)

where the p_k are the diagonal terms of $\hat{ρ}$ ⁠, each of which represents the probability of the system being in the kth energy eigenlevel, ε_k is the energy eigenvalue of the kth eigenlevel, p is the vector ${p_{k}}, 〈 \cdot 〉$ is the expectation value of the property given $\hat{ρ}$ ⁠, and τ is the relaxation time. The equation of motion, which generally consists of an infinite number of ordinary differential equations, can be solved using the concept of hypoequilibrium state and/or the density of states method, both of which were developed by Li and von Spakovsky [16]. A brief introduction to the concept of hypoequilibrium state is given in Sec. 2.3.

Nonequilibrium Evolution: Kinetics and Dynamics.

In general, the equation of motion for a system with a given group of conservation laws takes the following form:

\frac{d p_{k}}{d t} = \frac{1}{τ (p)} D_{k} (p)

(4)

where D_k and τ can be a function of p. Defining a new parameter

\tilde{τ}

⁠, such that

d \tilde{τ} = \frac{1}{τ (p (t))} d t, or \tilde{τ} = \int_{path} \frac{1}{τ (p (t'))} d t' = \tilde{τ} (t)

(5)

where

\tilde{τ}

is called the dimensionless time, the independent variable for the equation of motion can be changed so that

\frac{d p_{k}}{d \tilde{τ}} = D_{k} (p)

(6)

The solution for this equation is written as

p_{k} = p_{k} (\tilde{τ})

(7)

No matter how the relaxation time τ depends on the real time t and the state, the equation of motion can always be transformed to Eq. (6) with the parameter change defined by Eq. (5). Furthermore, the evolution of system state follows the same function (Eq. (7)) in $\tilde{τ}$ ⁠. Physically, this means that it follows the same thermodynamic trajectory in state space.

By doing the preceding, the kinetics of the system evolution is found via Eqs. (6) and (7) and results in a unique trajectory in state space based on the parameter $\tilde{τ}$ ⁠. This trajectory gives the intermediate states through which the system passes during state relaxation and the sequence in which they occur. The dynamics is found via Eq. (4) and the functional dependence $τ = τ (p)$ and results in the same trajectory in state space but this time based on the real time t. In this paper, a constant relaxation time τ calculated via kinetic theory is used (see Sec. 3.2). A constant τ allows a discussion on the kinetics of the system state evolution, which can be regarded as a thermodynamic property of the system's nonequilibrium evolution. Generally, the relaxation time τ is determined based on the microscopic interactions or phenomena causing the relaxation [17]. However, no matter what kind of microscopic interaction causes the relaxation, a system's nonequilibrium state evolution follows the same trajectory, one decided by its eigenstructure and the geometry of its system state space. In this way, the thermodynamic features of relaxation are revealed without going into the detailed microscopic mechanics of the system.

Nonequilibrium State and State Evolution Description: Hypoequilibrium.

The thermodynamic features of the nonequilibrium relaxation process generated by the SEAQT framework have a number of useful characteristics, which allow a unique and complete description of each nonequilibrium state via a set of extensive and intensive properties (as opposed to via a distribution function). Moreover, the definition of each nonequilibrium intensive property is fundamental relative to the laws of thermodynamics rather than phenomenological. Such a description is based on the key concept of hypoequilibrium state developed by Li and von Spakovsky [16,17], which is a direct result of the SEA principle. The physical meaning of hypoequilibrium state is that for a system in a nonequilibrium state, if a subset of its energy eigenlevels are in mutual equilibrium (i.e., the probability distribution for these levels is a Maxwellian distribution) that subset of levels will remain in mutual equilibrium throughout the entire state evolution of the system (i.e., the probability distribution of that subset remains Maxwellian). Thus, a process of relaxation can be described as follows. For any initial state (a distribution of probabilities among the energy eigenlevels), the eigenlevels can always be regrouped into many subsets such that in every subset the probability distribution is a Maxwellian distribution with a given intensive property (e.g., temperature) [16,20]. However, the different subsets are not necessarily in mutual equilibrium with each other. During the system relaxation process, the energy eigenlevels in one subset evolve together to reach mutual equilibrium with the other subsets. The concept of hypoequilibrium state is, thus, well defined for any state of the system and is fundamental and, therefore, much more general than the local equilibrium assumption. The latter discretizes a global system into many small local systems, each of which is assumed at a phenomenological level to approximately be in a state of equilibrium. In contrast, the hypoequilibrium concept permits each local system to be represented by a nonequilibrium state (even on far-from-equilibrium) for which the probability distribution is a non-Maxwellian distribution. In this section, the mathematical representation of hypoequilibrium is given using a system whose state can be characterized by energy and particle number and a Maxwellian distribution with the intensive properties of temperature and chemical potential. For proofs and a more detailed discussion on the hypoequilibrium state concept, the reader is referred to Refs. [16,17], and [20].

For a given system represented by an energy eigenlevel set

Ω = {(ε_{k}, N_{k})}

⁠, where each energy eigenlevel is represented by a pair of energy and particle number eigenvalues, the system can be divided into M subsystems

Ω_{i} = {(ε_{k}^{i}, N_{k}^{i})}, Ω_{i} \cap Ω_{j} = \emptyset

and

Ω = \cup Ω_{i}

⁠. For a system with E and N conserved (or only E conserved), if the probability distribution in each subsystem yields to a grand canonical distribution (or canonical distribution), the system is designated as being in an Mth-order hypoequilibrium state. A hypoequilibrium state of order 1 corresponds to a state in stable equilibrium. For a system with E and N conserved, the probability distribution of the Mth-order hypoequilibrium state takes the following form:

\forall i = 1, 2, ..., M, p_{k}^{i} = \frac{p^{i}}{Ξ^{i} (β^{i}, γ^{i})} e^{- β^{i} ε_{k}^{i} - γ^{i} N_{k}^{i}}, k = 1, 2, ..., w_{i}

(8)

where βⁱ and γⁱ are the parameters, pⁱ is the total probability in subsystem i, and

Ξ^{i} (β^{i}, γ^{i})

is the grand partition function of the subsystem with parameters βⁱ and γⁱ. The grand partition function is written as

Ξ^{i} (β^{i}, γ^{i}) = \sum_{k = 1}^{w_{i}} e^{- β^{i} ε_{k}^{i} - γ^{i} N_{k}^{i}}

(9)

For a given Mth-order hypoequilibrium state, the intensive properties of the subsystems can be represented by βⁱ and γⁱ, or equivalently using the temperature Tⁱ and chemical potential μⁱ defined by

T^{i} = \frac{1}{k_{b} β^{i}}, μ^{i} = γ^{i} T^{i}

(10)

A Mth-order hypoequilibrium state can then be represented by a division

{Ω_{i} = (ε_{k}^{i}, N_{k}^{i}), i = 1, ..., M}

of the system and a corresponding triplet set

{(p^{i}, β^{i}, γ^{i}), i = 1, ..., M}

⁠. The intensive property set

{(T^{i}, μ^{i}), i = 1, ..., M}

is a generalization of the definition of intensive property at stable equilibrium

(T^{eq}, μ^{eq})

⁠, which is also a first-order hypoequilibrium state. Li and von Spakovsky proved in Ref. [16] that if a system begins in an Mth-order hypoequilibrium state, it will remain in an Mth-order hypoequilibrium state throughout the state evolution as long as the same subsystem division is maintained. Thus, the time evolution of the distribution takes the following form:

\begin{matrix} \forall i = 1, 2, ..., M, p_{k}^{i} (t) = \frac{p^{i} (t)}{Ξ^{i} (β^{i} (t), γ^{i} (t))} e^{- β^{i} (t) ε_{k}^{i} - γ^{i} (t) N_{k}^{i}} \\ , k = 1, 2, ..., w_{i} \end{matrix}

(11)

The intensive property set ${(T^{i} (t), μ^{i} (t)), i = 1, ..., M}$ is also well defined throughout the entire evolution, and the system evolution can also be represented by the evolution of the triplet set ${(p^{i} (t), β^{i} (t), γ^{i} (t)), i = 1, ..., M}$ for a total of only 3 M variables.

SEAQT for Interacting Systems

Heat and Mass Interactions.

The equation of motion is designed to study the nonequilibrium relaxation process of an isolated system. However, since interacting systems can be viewed as a composite isolated nonequilibrium system of subsystems whose interactions cause the relaxation process of the composite, the SEAQT equation of motion calculate can be used to determine the evolution of these interacting systems as well. Furthermore, if the states of the subsystems of the composite are themselves not represented by a grand canonical (or canonical) distribution (i.e., they are in nonequilibrium states), higher order hypoequilibrium states [16,20] can be applied and used to study the composite system and its subsystem interactions.

The SEAQT framework, thus, explains these interactions on the basis of thermodynamics. This contrasts with other nonequilibrium approaches, which attempt to predict the nonequilibrium relaxation process on the basis of mechanical microscopic interactions. It also contrasts with traditional thermodynamics, which necessarily limits its description to quasi-equilibrium processes and the near-equilibrium realm. The SEAQT framework, on the other hand, is able to not only describe these processes and this realm but also nonquasi-equilibrium processes (e.g., see Ref. [18]) and the far-from-equilibrium realm as well [14–20].

Now, consider two interacting systems (systems a and b). The composite system is not in stable equilibrium, and its state is represented by the probability distributions

{p_{i}^{a}}

in system a and

{p_{i}^{b}}

in system b together. There are four conservation laws: probability normalization of system a, probability normalization of system b, energy conservation of the composite system, and particle number conservation of the composite system. Based on SEA and the conservation laws,

{{\hat{I}}^{a}, {\hat{I}}^{b}, \hat{H}, \hat{N}}

serve as the generators of the motion and the equation of motion for system a takes the following form [20,23]:

\frac{d p_{k}^{a}}{d t} = \frac{1}{τ} \frac{| \begin{matrix} - p_{k}^{a} \ln p_{k}^{a} & p_{k}^{a} & 0 & ε_{k}^{a} p_{k}^{a} & N_{k}^{a} p_{k}^{a} \\ {〈 s 〉}^{a} & 1 & 0 & {〈 e 〉}^{a} & {〈 N 〉}^{a} \\ {〈 s 〉}^{b} & 0 & 1 & {〈 e 〉}^{b} & {〈 N 〉}^{b} \\ 〈 e s 〉 & {〈 e 〉}^{a} & {〈 e 〉}^{b} & 〈 e^{2} 〉 & 〈 e N 〉 \\ 〈 N s 〉 & {〈 N 〉}^{a} & {〈 N 〉}^{b} & 〈 e N 〉 & 〈 N^{2} 〉 \end{matrix} |}{| \begin{matrix} 1 & 0 & {〈 e 〉}^{a} & {〈 N 〉}^{a} \\ 0 & 1 & {〈 e 〉}^{b} & {〈 N 〉}^{b} \\ {〈 e 〉}^{a} & {〈 e 〉}^{b} & 〈 e^{2} 〉 & 〈 e N 〉 \\ {〈 N 〉}^{a} & {〈 N 〉}^{b} & 〈 e N 〉 & 〈 N^{2} 〉 \end{matrix} |}

(12)

where

{〈 \cdot 〉}^{a (b)}

is the expectation value of a property in system a(b), and

〈 \cdot 〉 = {〈 \cdot 〉}^{a} + {〈 \cdot 〉}^{b}

is the expectation value of the total property of the composite system. The numerator of the fraction to the right of the equals can be expanded by the elements of the first row and their cofactors such that

\det = - p_{k}^{a} \ln p_{k}^{a} C_{1} + p_{k}^{a} C_{2}^{a} + ε_{k}^{a} p_{k}^{a} C_{3} + N_{k}^{a} p_{k} C_{4}

(13)

where C₁, C₂, C₃, and C₄ are the cofactors of the first line of the determinant. By defining

\frac{C_{2}^{a}}{C_{1}} = - α^{a}, \frac{C_{3}}{C_{1}} = - β, \frac{C_{4}}{C_{1}} = - γ

(14)

the equation of motion transforms to

\frac{d p_{k}^{a}}{d t} = \frac{1}{τ} (- p_{k}^{a} \ln p_{k}^{a} - p_{k}^{a} α^{a} - ε_{k}^{a} p_{k}^{a} β - N_{k}^{a} p_{k}^{a} γ)

(15)

If the two systems are in stable equilibrium states, the composite system is in a second-order hypoequilibrium state. Li and von Spakovsky [20] have proven that the conclusions in Sec. 2.3 relative to the hypoequilibrium state still hold for the composite system. Thus, the time evolution of the two systems takes the following form:

System a (b) : p_{k}^{a (b)} (t) = \frac{e^{- γ^{a (b)} (t) N_{k}^{a (b)} - β^{a (b)} (t) ε_{k}^{a (b)}}}{Ξ^{a (b)} (β^{a (b)} (t), γ^{a (b)} (t))}

(16)

The time evolution of

β^{a} (t), β^{b} (t), γ^{a} (t)

⁠, and

γ^{b} (t)

can be solved from the equation of motion for intensive properties [20], such that

system a (b) : \frac{d β^{a (b)} (t)}{d t} = \frac{1}{τ} (β^{a (b)} (t) - β (t))

(17)

system a (b) : \frac{d γ^{a (b)} (t)}{d t} = \frac{1}{τ} (γ^{a (b)} (t) - γ (t))

(18)

where β and γ are given by Eq. (14). For the general form where two systems are in nonequilibrium states, the reader is referred to the detailed discussion given in Ref. [20].

Relaxation Time τ.

Depending on which type of interaction is of interest, the relaxation time τ can be estimated using fundamental transport information such as a diffusion coefficient, a chemical reaction rate constant [16,19], or a viscosity. The self-diffusion coefficient is used here. If we assumed that diffusion dominates the relaxation of the system, the relation between the system relaxation time and its self-diffusion coefficient D_self as derived in Ref. [16] is given by

τ = \frac{{(δ x)}^{2}}{2 D_{self}}

(19)

where δx is the distance between the centers of the interacting systems (approximated by the dimension of system a here) and

D_{self} = D_{self} (P, T)

is a function of the pressure P and the temperature T. The dependence of the diffusion coefficient on temperature for gases is acquired by using the Chapman–Enskog theory. For the binary diffusion of A and B, the diffusion coefficient takes the form

D_{A B} = \frac{0.0018583 T^{3 / 2} \sqrt{1 / M_{A} + 1 / M_{B}}}{P σ_{A B}^{2} Ω}

(20)

where $P [atm]$ is the pressure, $T [K]$ is the temperature, M_A and M_B are the standard atomic weight, $σ_{A B} [Å]$ is the average collision diameter, and Ω is a temperature-dependent collision integral. In the study below of the self-diffusion of Helium-3, the parameters of A and B are both chosen to be that of the Helium-3 molecule, while $σ_{A B} [Å]$ and Ω are both taken from tabulated data [24]. For simplicity, the self-diffusion coefficient for Helium-3 is used to get the order of magnitude for a constant τ. For Helium-3 at 1 atm and 300 K, $τ \sim 10^{- 5}$ s for $δ x = 1$ mm. In order to study the indistinguishable effect of Helium-3 and Helium-4 on the thermodynamic trajectories predicted minus the influence of the isotope mass effect on relaxation time (see Eq. (20)), dimensionless time is used in the results.

Many Particle System.

For a many particle system, the occupation number representation is used, and the state space is a Fock space, which is the sum of N-particle state spaces [25]. The N-particle basis state is

| n_{ν_{1}}, n_{ν_{2}}, n_{ν_{3}}, \dots 〉, \sum_{k} n_{ν_{k}} = N

⁠. The space spanned by the occupation number basis is the Fock space

F = \oplus_{N = 0}^{\infty} F_{N}

⁠, where

F_{N} = span {| n_{ν_{1}}, n_{ν_{2}}, n_{ν_{3}}, \dots 〉 | \sum_{k} n_{ν_{k}} = N}

⁠. For independent particles, the occupation number basis is also the energy eigenlevels of the Hamiltonian. The occupation number

n_{ν_{k}}

is the particle number distributed in single-particle energy eigenlevels ν_k with values of

fermions : n_{ν_{k}} = 0, 1

(21)

bosons : n_{ν_{k}} = 0, 1, 2, \dots

(22)

In order to get system properties, the grand partition function is calculated, which is defined by

Ξ (β, γ) = \sum_{N} e^{- γ N} \sum_{S_{N}} e^{- β E_{S_{N}}}

(23)

where

E_{S_{N}}

is the N-particle energy eigenvalue of the S_Nth eigenlevel. The sum over S_N is a summation over N-particle energy eigenlevels. The natural logarithm of the partition function in terms of the single-particle energy eigenlevels ε_k then takes the following form:

\ln Ξ (β, γ) = \pm \sum_{k} \ln (1 \pm e^{- γ - β ε_{k}})

(24)

where fermions take the plus sign and bosons the minus sign.

The single-particle energy eigenlevels are chosen to be those for free particles in a box with Helium-3 as fermions and Helium-4 as bosons. Thus,

ε (n_{x}, n_{y}, n_{z}) = \frac{ℏ^{2}}{8 m} ({(\frac{n_{x}}{L_{x}})}^{2} + {(\frac{n_{y}}{L_{y}})}^{2} + {(\frac{n_{z}}{L_{z}})}^{2})

(25)

where n_x, n_y, and n_z are the quantum numbers for the translational energy eigenlevels, L_x, L_y, and L_z are the three dimensions of the box,

ℏ

is the reduced Plank constant, and m is the mass of a particle. The volume is given by

V = L_{x} L_{y} L_{z}

⁠. If the system yields to the two conditions

e^{- γ} < 1, β ≪ \frac{1}{Δ ε} = \frac{8 m L^{2}}{ℏ^{2}}

(26)

which mean that the concentration is not too high (i.e., the weak degenerate condition) and the temperature is much higher than the characteristic temperature for translation (i.e., the quasi-continuous condition). The natural logarithm of the partition function can then be calculated using

\ln Ξ (β, γ) = \mp \frac{V}{λ_{T}^{3}} (2 s + 1) L i_{5 / 2} (\mp e^{- γ})

(27)

where fermions take the minus sign and bosons the plus sign, s is the spin of a particle,

(2 s + 1)

is the degeneracy of the spin,

L i_{s} (z)

is the polylogarithm function defined by

L i_{s} (z) = \sum_{k}^{\infty} \frac{z^{k}}{k^{s}}

(28)

and λ_T, a function of β, is the de Broglie wavelength of the thermal energy defined by

λ_{T} = \frac{h}{\sqrt{2 π m k_{b} T}} = \frac{h}{\sqrt{2 π m / β}}

(29)

The properties for system a are then given by

{〈 N 〉}^{a} = \mp \frac{V}{λ_{T}^{3}} (2 s + 1) L i_{3 / 2} (\mp e^{- γ^{a}})

(30)

{〈 e 〉}^{a} = \mp \frac{V}{λ_{T}^{3}} \frac{3}{2 β^{a}} (2 s + 1) L i_{5 / 2} (\mp e^{- γ^{a}})

(31)

{〈 s 〉}^{a} = k_{b} (\frac{5}{2} \ln Ξ^{a} + γ^{a} {〈 N 〉}^{a})

(32)

{〈 e^{2} 〉}^{a} = {({〈 E 〉}^{a})}^{2} + \frac{15}{4 {(β^{a})}^{2}} \ln Ξ^{a}

(33)

{〈 N^{2} 〉}^{a} = {({〈 N 〉}^{a})}^{2} \mp \frac{V}{λ_{T}^{3}} (2 s + 1) L i_{1 / 2} (\mp e^{- γ^{a}})

(34)

{〈 e N 〉}^{a} = {〈 E 〉}^{a} {〈 N 〉}^{a} \mp \frac{V}{λ_{T}^{3}} \frac{3}{2 β^{a}} (2 s + 1) L i_{3 / 2} (\mp e^{- γ^{a}})

(35)

{〈 s N 〉}^{a} = k_{b} (β^{a} {〈 e N 〉}^{a} + γ^{a} {〈 N^{2} 〉}^{a} + \ln Ξ^{a} {〈 N 〉}^{a})

(36)

{〈 e s 〉}^{a} = k_{b} (β^{a} {〈 e^{2} 〉}^{a} + γ^{a} {〈 e N 〉}^{a} + \ln Ξ^{a} {〈 e 〉}^{a})

(37)

where fermions take the minus sign and bosons the plus sign. The properties for system b have the same form. These properties are used to calculate the β and γ in Eq. . From Eq. , one can arrive at the expression for the particle number concentration, namely,

n^{a (b)} = \mp \frac{(2 s + 1) L i_{3 / 2} (\mp e^{- γ^{a (b)}})}{λ_{T}^{3}}

(38)

where again fermions take the minus sign and bosons the plus sign. It should be noted that the restrictions on concentration and temperature (Eq. (26)) employed here are only used to obtain a simple expression for the partition function (Eq. (27)) and do not in anyway limit the generality of the theory nor the thermodynamic conclusions described in Secs. 2 and 3.1. When condition (26) does not hold (e.g., in the case of Bose–Einstein condensation), one can evaluate the partition function using Eq. (24) and apply the SEAQT framework in the same way illustrated previously. Applications of this framework for more general conditions will be presented in future studies.

Results and Discussion

The SEAQT model is applied to study the size effects on heat and mass diffusion. The discussion of size effects focuses on their influence on the thermodynamic properties or kinetics of the relaxation process since the results are shown in dimensionless time. There are two ways for quantum mechanics to influence the thermodynamic properties of the system. One is from the fact that energy eigenlevels are discrete at temperatures close to or lower than the characteristic temperature for a given energy mode (e.g., that of translation, rotation, vibration, etc.). The other is from the indistinguishable feature of a particle, where statistical correlations are important. The results presented here are limited to the realm where Eq. (26) holds, i.e., where the temperature is significantly higher than the characteristic temperature for translation (condition #2 of Eq. (26)) and the chemical potential is not too low (condition #1 of Eq. (26)) so that the system is in the weak degenerate realm. For strong degenerate effects, such as with Bose–Einstein condensation (i.e., $e^{γ} ≪ 1$ or $e^{- γ} ≫ 1$ ⁠), Eq. (27) needs an additional term for the energy ground state, but this is left for a future paper. Nonetheless, the indistinguishable feature of the particles still influences the nonequilibrium relaxation process in this realm when $e^{γ}$ is not too large (⁠ $e^{γ} < 10^{3}$ ⁠) or $e^{- γ}$ is not too small (⁠ $e^{- γ} > 10^{- 3}$ ⁠). Furthermore, since the temperature is much higher than the characteristic translational temperature, the quantum effect from discrete energy eigenlevels is not studied here and is also left for a future paper. For the discrete case, the summation in Eq. (24) could not be calculated by the polylogarithm function.

The particle number evolutions and temperature evolutions for bosons and fermions are shown in Figs. 1–4. The temperatures of system a and system b are chosen to be 300 K and 1000 K, respectively. In Figs. 1 and 2, the $e^{γ}$ of system a and system b are chosen to be 1.2 and 1.5, respectively. In Figs. 3 and 4, the $e^{γ}$ of system a and system b are chosen to be 1.2 × 10⁵ and 1.5 × 10⁵, respectively. The solid line is for the interacting boson systems, and the dashed line is for the interacting fermion systems. In order to make a comparison, in all the four figures (Figs. 1–4), fermions and bosons are assumed to have the same particle mass $m_{{He}_{3}}$ and the same box size of 10^–3 m^–3. This comparison also necessitates removing the degenerate influence, which requires that the particle number, Eq. (30), for system a and the particle number for system b be divided by their respective degeneracies, which for a single particle with spin s is $2 s + 1$ ⁠. From these four figures, one can observe that the difference in the behavior of fermions and bosons decreases as $e^{γ}$ increases which in turn means that the concentration (Eq. (38)) decreases. At the low concentration limit, $e^{γ} ≫ 1$ ⁠, both fermions and bosons go back to being classical particles.

Fig. 1

Temperature evolutions for the boson systems (solid line) and that for the fermion systems (dashed line). The concentration is high, because eγ∼10.

View large Download slide

Temperature evolutions for the boson systems (solid line) and that for the fermion systems (dashed line). The concentration is high, because $e^{γ} \sim 10$ ⁠.

Fig. 2

Particle number evolutions for the boson systems (solid line) and that for the fermion systems (dashed line). The concentration is high, because eγ∼10.

View large Download slide

Particle number evolutions for the boson systems (solid line) and that for the fermion systems (dashed line). The concentration is high, because $e^{γ} \sim 10$ ⁠.

Fig. 3

View large Download slide

Temperature evolutions for the boson systems (solid line) and that for the fermion systems (dashed line). The solid line and the dashed line converge. The concentration is low, because $e^{γ} ≫ 1$ ⁠.

Fig. 4

View large Download slide

Particle number evolutions for the boson systems (solid line) and that for fermion systems (dashed line). The solid line and the dashed line converge. The concentration is low, because $e^{γ} ≫ 1$ ⁠.

In Figs. 5 and 6, four cases for the boson systems are modeled in order to study the effect due to concentration changes. The temperatures of the two systems are those for Figs. 1–4. $(e^{γ^{a}}, e^{γ^{a}})$ are chosen to be $(1.2, 1.5)$ ⁠, (12, 15), (1200, 1500), and $(1.2 \times 10^{5}, 1.5 \times 10^{5})$ ⁠, where the difference $γ^{a} - γ^{b}$ remains the same. It can be observed that the concentration effect influences the thermodynamic trajectory of the nonequilibrium evolutions via changes in γ (chemical potential) or equivalently via changes in concentration (Eq. (38)). Systems with the same particle number but different concentrations follow this type of behavior.

Fig. 5

View large Download slide

Temperature evolutions for the boson systems for the four cases. γ^a and γ^b increase and, thus, concentration decreases from the bottom curve to the top one for system a and from the top curve to the bottom one for system b. $γ^{a} - γ^{b}$ is kept constant. The curves for the two low concentration cases converge.

Fig. 6

View large Download slide

Normalized particle number evolutions for the boson systems for the four cases. γ^a and γ^b increase and, thus, concentration decreases from the bottom curve to the top one for system a and from the top curve to the bottom one for system b. $γ^{a} - γ^{b}$ is kept constant. The curves for the two low concentration cases converge.

The size or volume effect can be studied by looking at systems with the same concentration but different volumes. Figures 7 and 8 show the particle number evolutions and temperature evolutions for boson system interactions for the different volumes but fixed γ^a and γ^b (or equivalently, system concentrations) and other parameters. The particle number curves are normalized by the total particle number at any instant of time. The $e^{γ}$ of systems a and b are chosen to be 1.2 and 1.5, respectively, and their temperatures 300 K and 1000 K. The volumes of the systems are chosen to be $10^{- 3}$ m³, $10^{- 9}$ m³, and $10^{- 18}$ m³. It can be observed that there is no difference in normalized trajectory. This can be proven by substituting Eqs. (30)–(37) into Eq. (12) to acquire the explicit expression for β and γ. Actually, under the conditions of Eq. (26), β and γ do not depend on V explicitly. Thus, the trajectory is not influenced by volume. However, for the case when the conditions of Eq. (26) do not hold, the natural logarithm of the grand partition function no longer depends linearly on volume as it does in Eq. (27). In this case, the size effect may appear and the summation of Eq. (24) cannot be calculated by polylogarithm function. Moreover, the size effect can also appear when the study is conducted in real time. However, this volume effect can only influence the speed of the system going through the trajectory but cannot change the trajectory or intermediate states of system state evolution.

Fig. 7

View large Download slide

Temperature evolutions for the boson systems for different volumes and the same $(γ^{a}, γ^{b})$ ⁠. The curves of the three cases converge.

Fig. 8

View large Download slide

Normalized particle number evolution evolutions for boson systems for different volumes and the same $(γ^{a}, γ^{b})$ ⁠. The curves of the three cases converge.

Conclusions

The SEAQT framework is a powerful and novel approach for studying nonequilibrium phenomena, even those far-from-equilibrium, from the atomistic level up to the macroscopic. With this framework, heat and mass diffusion, even in realms far-from-equilibrium, can be modeled. Furthermore, with the concepts of hypoequilibrium state and nonequilibrium intensive property, a description of nonequilibrium states comparable to that at stable equilibrium is now possible. To illustrate this framework and these concepts, the kinetics of the nonequilibrium evolutions of heat and mass diffusion were predicted for interacting boson as well as fermion systems, and their state evolution trajectories in thermodynamic state space are determined. Emphasis is placed on the thermodynamic features of the trajectories with each revealing through which intermediate thermodynamic states the system evolves and by which sequence these states occur.

In addition, two kinds of effects influencing nonequilibrium trajectories have been studied. The concentration effect is shown to influence the state evolution trajectory in the weak degenerate realm, which comes from the indistinguishable feature of the particles. In contrast, size effect shows no influence on these trajectories in the high temperature limit where the discrete energy eigenlevels are treated as continuous. This is the case even though these changes influence system properties via the grand partition function.

Acknowledgment

The funding for this research was provided by the U.S. Office of Naval Research under ONR Grant No. N00014-11-1-0266.

Nomenclature

$\hat{D}$ =: dissipation term in the SEAQT equation of motion
D_self =: phenomenological self-diffusion coefficient
e =: energy of the system
$\hat{H}$ =: Hamiltonian
$L i_{S} (z)$ =: polylogarithm function
M =: order of the hypoequilibrium state, i.e., the number of subsystems
$n^{a (b)}$ =: particle number concentration in the system a(b)
N_k =: particle number of the kth eigenlevel
p =: ${p_{k}}$
p_k =: probability distribution among energy eigenlevels
$p_{k}^{a (b)}$ =: probability distribution among energy eigenlevels in the system a(b)
s =: entropy of the system
Tⁱ =: temperature of the ith subsystem
V =: volume of the system
$X_{θ}^{H}$ =: reversible term in the Ginzburg–Landau equation
$Y_{θ}^{H}$ =: reversible term in the Ginzburg–Landau equation
$Y_{θ (t)}^{H}$ =: irreversible term in the Ginzburg–Landau equation
$α^{a (b)}, β, γ$ =: intermediate parameters of the equation of motion in the system a(b), defined by Eq. (14)
$β^{a (b)}$ =: intensive properties corresponding to the energy in the system a(b)
βⁱ =: intensive properties corresponding to the energy in the ith subsystem
$γ^{a (b)}$ =: intensive properties corresponding to the particle number in the system a(b)
γⁱ =: intensive properties corresponding to the particle number in the ith subsystem
δx =: dimension of system a
ε_k =: energy of the kth eigenlevel
$ε_{k}^{a (b)}$ =: energy of the kth eigenlevel in the system a(b)
$θ (t)$ =: state evolution trajectory in the Ginzburg–Landau equation
λ_T =: de Broglie wavelength
μⁱ =: chemical potential of the ith subsystem
$Ξ^{i}$ =: partition function in the ith subsystem and a function of $(β^{i}, γ^{i})$
$\hat{ρ}$ =: density operator
τ =: relaxation time in the SEAQT equation of motion
$\tilde{τ}$ =: dimensionless time
Ω =: the set of the eigenlevels in the system
Ω_i =: the set of the eigenlevels in the ith subsystem
$〈 \cdot 〉$ =: expectation value
${〈 \cdot 〉}^{a (b)}$ =: expectation value in the system a(b)

References

1.

De Groot

,

S. R.

, and

Mazur

,

P.

,

2013

,

Non-Equilibrium Thermodynamics

,

Courier

,

New York

.

2.

Zanchini

,

E.

, and

Beretta

,

G. P.

,

2014

, “

Recent Progress in the Definition of Thermodynamic Entropy

,”

Entropy

,

16

(

3

), pp.

1547

–

1570

.

Google Scholar

Crossref

3.

Grmela

,

M.

, and

Öttinger

,

H. C.

,

1997

, “

Dynamics and Thermodynamics of Complex Fluids—I: Development of a General Formalism

,”

Phys. Rev. E

,

56

(

6

), pp.

6620

–

6632

.

Google Scholar

Crossref

4.

Öttinger

,

H. C.

, and

Grmela

,

M.

,

1997

, “

Dynamics and Thermodynamics of Complex Fluids—II: Illustrations of a General Formalism

,”

Phys. Rev. E

,

56

(

6

), pp.

6633

–

6655

.

Google Scholar

Crossref

5.

Chen

,

S.

, and

Doolen

,

G. D.

,

1998

, “

Lattice Boltzmann Method for Fluid Flows

,”

Ann. Rev. Fluid Mech.

,

30

(

1

), pp.

329

–

364

.

Google Scholar

Crossref

6.

Rapaport

,

D. C.

,

2004

,

The Art of Molecular Dynamics Simulation

,

Cambridge University Press

,

Cambridge, UK

.

7.

Hoffmann

,

K. H.

,

2001

, “

Quantum Thermodynamics

,”

Ann. Phys.

,

10

(

1–2

), pp.

79

–

88

.

Google Scholar

Crossref

8.

Beretta

,

G. P.

,

Gyftopoulos

,

E. P.

,

Park

,

J. L.

, and

Hatsopoulos

,

G. N.

,

1984

, “

Quantum Thermodynamics. A New Equation of Motion for a Single Constituent of Matter

,”

Nuovo Cimento B

,

82

(

2

), pp.

169

–

191

.

Google Scholar

Crossref

9.

Beretta

,

G. P.

,

Gyftopoulos

,

E. P.

, and

Park

,

J. L.

,

1985

, “

Quantum Thermodynamics. A New Equation of Motion for a General Quantum System

,”

Nuovo Cimento B

,

87

(

1

), pp.

77

–

97

.

Google Scholar

Crossref

10.

Beretta

,

G.

,

2006

, “

Nonlinear Model Dynamics for Closed-System, Constrained, Maximal-Entropy-Generation Relaxation by Energy Redistribution

,”

Phys. Rev. E

,

73

(

2

), p.

026113

.

Google Scholar

Crossref

11.

Beretta

,

G. P.

,

2010

, “

Maximum Entropy Production Rate in Quantum Thermodynamics

,”

J. Phys.: Conf. Ser.

,

237

(1), p.

012004

.

Google Scholar

Crossref

12.

Beretta

,

G. P.

,

2014

, “

Steepest Entropy Ascent Model for Far-Nonequilibrium Thermodynamics: Unified Implementation of the Maximum Entropy Production Principle

,”

Phys. Rev. E

,

90

(

4

), p.

042113

.

Google Scholar

Crossref

13.

von Spakovsky

,

M.

, and

Gemmer

,

J.

,

2014

, “

Some Trends in Quantum Thermodynamics

,”

Entropy

,

16

(

6

), pp.

3434

–

3470

.

Google Scholar

Crossref

14.

Li

,

G.

,

Al-Abbasi

,

O.

, and

von Spakovsky

,

M. R.

,

2014

, “

Atomistic-Level Non-Equilibrium Model for Chemically Reactive Systems Based on Steepest-Entropy-Ascent Quantum Thermodynamics

,”

J. Phys.: Conf. Ser.

,

538

(1), p.

012013

.

Google Scholar

Crossref

15.

Cano-Andrade

,

S.

,

Beretta

,

G. P.

, and

von Spakovsky

,

M. R.

,

2015

, “

Steepest-Entropy-Ascent Quantum Thermodynamic Modeling of Decoherence in Two Different Microscopic Composite Systems

,”

Phys. Rev. A

,

91

(

1

), p.

013848

.

Google Scholar

Crossref

16.

Li

,

G.

, and

von Spakovsky

,

M. R.

,

2016

, “

Steepest-Entropy-Ascent Quantum Thermodynamic Modeling of the Relaxation Process of Isolated Chemically Reactive Systems Using Density of States and the Concept of Hypoequilibrium State

,”

Phys. Rev. E

,

93

(

1

), p.

012137

.

Google Scholar

Crossref

PubMed

17.

Li

,

G.

, and

von Spakovsky

,

M. R.

,

2016

, “

Generalized Thermodynamic Relations for a System Experiencing Heat and Mass Diffusion in the Far-From-Equilibrium Realm Based on Steepest Entropy Ascent

,”

Phys. Rev. E

,

94

(

3

), p.

032117

.

Google Scholar

Crossref

PubMed

18.

Li

,

G.

, and

von Spakovsky

,

M. R.

,

2016

, “

Modeling the Nonequilibrium Effects in a Nonquasi-Equilibrium Thermodynamic Cycle Based on Steepest Entropy Ascent and an Isothermal-Isobaric Ensemble

,”

Energy

,

115

(

Pt. 1

), pp.

498

–

512

.

Google Scholar

Crossref

19.

Li

,

G.

, and

von Spakovsky

,

M. R.

,

2015

, “

Study of the Transient Behavior and Microstructure Degradation of a SOFC Cathode Using an Oxygen Reduction Model Based on Steepest-Entropy-Ascent Quantum Thermodynamics

,”

ASME

Paper No. IMECE2015-53726.

20.

Li

,

G.

, and

von Spakovsky

,

M. R.

,

2016

, “

Generalized Thermodynamic Definitions of Reservoir and Measurement Based on Interacting Nonequilibrium Systems in the Far-From-Equilibrium Realm and the Principle of Steepest Entropy Ascent

,” e-print arXiv: 1601.02703.

21.

Grmela

,

M.

,

2014

, “

Contact Geometry of Mesoscopic Thermodynamics and Dynamics

,”

Entropy

,

16

(

3

), pp.

1652

–

1686

.

Google Scholar

Crossref

22.

Montefusco

,

A.

,

Consonni

,

F.

, and

Beretta

,

G. P.

,

2015

, “

Essential Equivalence of the General Equation for the Nonequilibrium Reversible-Irreversible Coupling (Generic) and Steepest-Entropy-Ascent Models of Dissipation for Nonequilibrium Thermodynamics

,”

Phys. Rev. E

,

91

(

4

), p.

042138

.

Google Scholar

Crossref

23.

Beretta

,

G. P.

,

2009

, “

Nonlinear Quantum Evolution Equations to Model Irreversible Adiabatic Relaxation With Maximal Entropy Production and Other Nonunitary Processes

,”

Rep. Math. Phys.

,

64

(

1–2

), pp.

139

–

168

.

Google Scholar

Crossref

24.

Bird

,

R. B.

,

Stewart

,

W. E.

, and

Lightfoot

,

E. N.

,

2002

,

Transport Phenomena

,

Wiley

,

New York

.

25.

Bruus

,

H.

, and

Flensbergy

,

K.

,

2004

,

Many-Body Quantum Theory in Condensed Matter Physics: An Introduction

,

Oxford University Press

,

Oxford, UK

.

Study of Nonequilibrium Size and Concentration Effects on the Heat and Mass Diffusion of Indistinguishable Particles Using Steepest-Entropy-Ascent Quantum Thermodynamics

Introduction

SEAQT Equation of Motion

SEAQT Equation of Motion.

Nonequilibrium Evolution: Kinetics and Dynamics.

Nonequilibrium State and State Evolution Description: Hypoequilibrium.

SEAQT for Interacting Systems

Heat and Mass Interactions.

Relaxation Time τ.

Many Particle System.

Results and Discussion

Conclusions

Acknowledgment

Nomenclature

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Study of Nonequilibrium Size and Concentration Effects on the Heat and Mass Diffusion of Indistinguishable Particles Using Steepest-Entropy-Ascent Quantum Thermodynamics

Introduction

SEAQT Equation of Motion

SEAQT Equation of Motion.

Nonequilibrium Evolution: Kinetics and Dynamics.

Nonequilibrium State and State Evolution Description: Hypoequilibrium.

SEAQT for Interacting Systems

Heat and Mass Interactions.

Relaxation Time τ.

Many Particle System.

Results and Discussion

Conclusions

Acknowledgment

Nomenclature

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only