GRASSMANN PHASE SPACE THEORY

for the

BEC/BCS CROSSOVER

in

COLD FERMIONIC ATOMIC GASES

B. J. Dalton

Centre for Quantum Technology Theory, Swinburne University of Technology,
Melbourne, Victoria 3122, Australia

##  1 Abstract

Grassmann Phase Space Theory (GSPT) is applied to the BEC/BCS crossover in
cold fermionic atomic gases and used to determine the evolution (over either
time or temperature) of the Quantum Correlation Functions (QCF) that specify:
(a) the positions of the spin up and spin down fermionic atoms in a single
Cooper pair and (b) the positions of the two spin up and two spin down
fermionic atoms in two Cooper pairs The first of these QCF is relevent to
describing the change in size of a Cooper pair, as the fermion-fermion
coupling constant is changed through the crossover from a small Cooper pair on
the BEC side to a large Cooper pair on the BCS side. The second of these QCF
is important for describing the correlations berween the positions of the
fermionic atoms in two Cooper pairs, which is expected to be small at the BEC
or BCS sides of the crossover, but is expected to be significant in the strong
interaction unitary regime, where the size of a Cooper pair is comparable to
the separation between Cooper pairs. In GPST the QCF are ultimately given via
the stochastic average of products of Grassman stochastic momentum fields, and
GPST shows that the stochastic average of the products of Grassman stochastic
momentum fields at a later time (or lower temperature) is related linearly to
the stochastic average of the products of Grassman stochastic momentum fields
at an earlier time (or higher temperature), and that the matrix elements
involved in the linear relations are all c-numbers. Expressions for these
matrix elements corresponding to a small time increment (or a small
temperature change) have been obtained analytically, providing the formulae
needed for numerical studies of the evolution that are planned for a future
publication. Various initial conditions are considered, including those for a
non-interacting fermionic gas at zero temperature and a high temperature gas
(where the effect of the interactions can be ignored in the initial state).
These would be relevant for studying the time evolution of the creation of a
Cooper pair when the interaction is switched on via Feshbach resonance
methods, or the evolution as the temperature is lowered for particular choices
of the fermion-fermion coupling constant, such as where the coupling constant
is very large - corresponding to the unitary regime. Full derivations of the
expressions have been presented in the Appendices.

##  2 Introduction

Research into ultracold atomic gases [1] has been a major activity since the
1990’s when Bose-Einstein condensation was achieved for bosonic atoms. Non-
interacting untrapped bosonic atoms at zero temperature form a BEC, with a
macroscopic occupancy of the lowest single particle energy state, since there
is no Pauli exclusion principle that forbids multiple occupancy. Since the
2000’s ultracold gases of fermionic atoms have been prepared, and here quite
different effects occur. Non-interacting untrapped fermionic atoms at zero
temperature form a Fermi gas, with each energy state only being occupied by
atoms with different spins due to the Pauli exclusion principle. Consequently,
energy states are filled up to the Fermi energy $E_{F}$, whose associated wave
number kF is proportional to the inverse of the average separation between the
atoms. In both cases the single particle states are plane waves with momentum
hk.

When interactions are taken into account there are features in common though –
for both types of atom the interaction potential only has a range b of
$O(10^{-9}m)$, whilst the typical separation between atoms in these dilute
gases is of $O(5x10^{-7}m)$ making it unlikely that more than two atoms
interact at once. The interaction between a pair of atoms determines the
scattering length $a$ – which is related to the collision cross section.
However, recent advances involving Feshbach resonances have made it possible
to modify the scattering length [2] (for example by varying a static magnetic
field that is coupled to the atoms so that an unbound state in an open
scattering channel coincides with a bound state in a closed channel) so that
its value ranges from negative to positive across the resonance – where (if
decay processes are ignored) it would become infinite. At this resonance (or
unitary) situation the strong interaction regime applies, with different
effects to the weak interaction regimes well away from resonance on either
side.

One of the interesting phenomena involving interacting ultracold fermionic
atomic gases is the BEC/BCS crossover [3], [4], [5], [6]. If the atoms are
essentially unconfined (as in a weak trapping potential) and at zero
temperature the weak interaction regime differs for bosons and fermions. For
the simplest bosonic case of spin $0$ atoms the picture is that each atom
behaves essentially independently, with the effect of the interatomic
potentials being reflected in a mean field based on the average density of all
the other atoms. The behaviour is treated using the Gross-Pitaevskii equation
[7]. However, for the simplest fermionic case of spin ${\frac{1}{2}}$ atoms
there is the possibility of two atoms with opposite spins pairing up and
behaving as a single entity. The behaviour can be described using BCS theory
[3], (also a mean field theory). Here though there are two extreme
possibilities depending on whether the two-body scattering length $a$ is
positive or negative. For $1/(k_{F}a)<<0$ superfluid BCS behaviour occurs
based on loosely bound Cooper pairs of atoms with opposite spins and momenta
(analogous to the Cooper pairs of electrons in metallic superconductors). For
$1/(k_{F}a)>>0$ a BEC forms based on tightly bound pairs of atoms with
opposite spins and momenta that constitute a spin $0$ bosonic molecule. These
regimes involve weakly interacting fermions. By using Fesbach resonance
methods, the BEC/BCS crossover from the BEC side to the BCS side can be
achieved experimentally, including the intermediate strong interaction (or
unitary) regime, where the scattering length becomes very large. The behaviour
of the quantum correlation functions (QCF) describing the positions of the
fermionic atoms in the Cooper pairs can be measured using methods such as
Bragg spectroscopy [8].

However, in terms of presenting a challenge for developing theories of
interacting ultracold fermionic atomic gases it is the strong interaction
regime (where the scattering length $|a|>>1/k_{F}$), that is currently of
great interest. Here the pair size is on the order of the interparticle
separation $1/k_{F}$ [3]. There are a number of methods that can be used to
treat interacting fermi gases. These include: (a) Diagrammatic methods [9],
[10], (b) Field theory techniques [11] (c) Quantum Monte Carlo methods [12],
[13], (d) Density functional methods [14] as well as (e) BEC theory [3] and
other (f) Variational approaches [15]. All of these methods have advantages
and disadvantages. BEC theory provides a simple over-view of the crossover,
but it is based on a variational approach in which the quantum state is
assumed to be pairs of opposite spins each described by the same single pair
wavefunction. This may not be a good approximation in the unitary regime.
Quantum Monte-Carlo methods are also based around assumed choices of
variational functions, and in general the fermion sign problem leads to
difficulties. To quote Ref. [13] ”A typical wavefunction or density matrix of
a many-fermion system has a complex nodal structure that prevents an easy
stochastic process from driving the system to its ground state or desired
excited state. This problem, originating from the Pauli exclusion principle of
fermions, is known as the fermion-sign problem, because crossing a node of the
wavefunction or density matrix changes its sign, and thus disallows a
straightforward weight interpretation of the wavefunction or density matrix.”.
Field theory methods involve considering re-normalisation issues, which makes
the approach rather complicated. Diagrammatic methods involve trying to decide
which Feynman diagrams are important and which can be neglected. Density
functional methods have been mainly used in atomic physics. Hence there is
scope for developing other methods for treating the many-body fermionic
systems that are involved in ultracold fermionic atomic gases.

Phase space methods (which were developed to treat topics in quantum optics)
can also be applied to treat problems in ultracold atomic gases. These methods
are described in a number of textbooks [16], [17], [18], [19]. Essentially the
quantum density operator is represented by a distribution function
(functional) which depends on phase space variables (fields) representing mode
(field) annihilation, creation operators, and the average value of physical
quantities are then given by phase space integrals (functional integrals)
involving the distribution function and with physical quantities replaced by
functions of the phase space variables. Several distinct distribution
functions have been introduced, based on whether the physical quantities
involve normally ordered, symmetrically ordered or anti-normally ordered
expressions involving the mode (field) annihilation, creation operaors. These
are the Glauber-Sudarshan $P$, Wigner $W$ and Husimi $Q$ functions
respectively. Often the phase space variables (field) for annihilation,
creation operators are related by complex conjugation, but treatments
involving double phase space distributions where this does not apply also
exist [20]. The evolution of the density operator is mapped onto a Fokker-
Planck equation (functional FPE) for the distribution function (functional).
It is then common to introduce stochastic phase space variables (fields) which
satisfy Ito (or Stratonovich) stochastic differential (field) equations that
are equivalent to the FPE (FFPE) and which involve stochastic Wiener
increments, such that the phase space integrals that give the average value
for physical quantities can then be replaced stochastic averages of functions
of the now-stochastic phase space variables (fields). One of the general
advantages of phase space theory is that it does not require any presumptions
regarding what the actual behaviour of the system will turn out to be. By
contrast, BEC theory is based on an assumed state vector in which all the
fermions are in Cooper pairs.

Since quantum optics involves the physics of photons - which are bosons - it
is not surprising that the application of phase space methods in ultracold
atomic gases has mainly been for bosonic atoms such as in the treatment of
Bose-Einsten Condensates. However, although less common, phase space methods
can also be applied to treat fermionic atoms - though the treatment must take
into account the fact that the mode (or field) annihilation, creation
operators obey anti-commutation rules rather than the commutation rules that
apply in the bosonic case. This means that the square of a fermionic mode
annihilation or creation operator is zero, so the question then arises as to
what type of phase space variable is best suited for the fermionic situation?
For bosons the natural choice is to use c-numbers for the phase space
variables or fields - since they commute on multiplication and their square is
non-zero. For fermions it may seem more natural to use numbers that have
properties such as the square being zero and having anti-commutation
properties in multiplication - just as for the mode or field operators they
represent. Such quantities are known in mathematics as Grassmann numbers, and
the first phase space theory for fermions was formulated by Cahill and Glauber
in terms of Grassmann phase space variables [21]. One of the first
applications followed shortly later [26] and was to momentum correlations
between fermions of opposite spins in a $1D$ trap.

There are however phase space theories for fermions which are based on
c-number quantities. One approach involves introducing Gaussian projectors
specified by c-number phase space variables [22], where for $n$ modes these
are $1+n(2n-1)$ variables that determine the covariance matrix involved in
defining the Gaussian projector. In the treatment by Rosales et al [23] the
density operator is represented by a $Q$ distribution function of the c-number
variables, where $Q$ is a positive quantity that can be interpreted as a
probability. The $Q$ function is defined on a phase space equivalent to the
space of real, anti-symmetric matrices. The approach has recently been applied
to Majorana fermions [24], and the Fermi-Hubbard Hamiltonian [25]. Here the
mode annihilation $\widehat{c}_{i}$ and creation $\widehat{c}_{i}^{{\dagger}}$
operators being replaced by Hermitian Marjorana operators
$(\widehat{c}_{i}+\widehat{c}_{i}^{{\dagger}})$ and
$-i(\widehat{c}_{i}-\widehat{c}_{i}^{{\dagger}})$.

Further developments of phase space theory for fermions based on Grassmann
variables are set out in following textbook [19] and articles [27], [28], [29]
The same general phase space theory features apply as described previously. It
is convenient to represent the quantum density operator for the system by a
Grassmann (un-normalised B) distribution function [26] in a phase space, but
again the phase space variables replacing the fermion annihilation and
creation operators are Grassmann variables. For $n$ modes there are $2n$ phase
space variables if a double space is used. Here the projector used in the
expression for the density operator is that involving fermion (Bargmann)
coherent states, which are given by similar expressions to the boson
(Bargmann) coherent states that are used in phase space theories for bosonic
systems. Experimentally measureable quantities for topics such as the BEC/BCS
crossover are given by quantum correlation functions (QCF), and these can be
related theoretically to Grassmann phase space integrals involving the
distribution function and with the fermion operators being replaced by phase
space variables. As previously, evolution equations (over time or temperature)
for the density operator lead to Fokker-Planck equations for the distribution
function. These in turn are replaced by Ito stochastic equations for
stochastic Grassmann phase space variables, and QCF now given by stochastic
averages of products of these stochastic variables. A key point is that the
stochastic averages for QCF at a later time (or temperature) can be shown to
be related linearly to such stochastic averages at an earlier time via
matrices that only involve c-numbers. Even though these matrices involve
stochastic quantities such as Wiener increments, they are not dependent on
Grassmann variables, which enables computations to be carried out without
having to represent Grassmann variables on the computer, and this makes
Grassmannn phase space theory numerically computable. The initial stochastic
averages of products are obtained from the initial density operator. It should
also be noted that although the distribution function (functionals) are only
quasi-probabilities (since they are Grassmann functions and not c-number
functions), the expressions for QCF and the mean values of physical quantities
has the same formal structure as in standard probability theory - with the
physical quantities given by Grassmann variables being weighted by the (quasi)
probability and then integrated over - albeit involving Grasssmann integrals
rather than c-number integrals.

The utility of the theory has been successfully tested [30] on a simple
fermion system that can be treated exactly by standard methods, namely
determining the coherence between two Cooper pair states functions for
interacting spin $\frac{1}{2}$ fermions. .

In the present paper, we now apply the Grassmann phase space theory to
treating the BEC/BCS crossover - a much more complex topic. The QCF that are
studied are those specifying the position of the spin up and spin down
fermions in a single Cooper pair, and those specifying the positions of the
two spin up and two spin down fermions in two Cooper pairs. We focus
particularly on treating how the size of a Cooper pair changes through the
crossover (BEC theory [3] indicates it changes smoothly), and on how the
correlations between two Cooper pairs change through the crossover (BEC theory
[3] indicates there is little correlation at the BEC or BCS ends, but strong
correlation should occur in the unitary regime). Both time evolution and
temperature evolution will be studied, the former from an initial state at
near zero temperature for non-interacting fermions with the fermions confined
within the Fermi sphere, the latter from an initial state at high temperature
where the fermion interactions can be ignored.

The present paper sets out the theory for the equations that determine the
above QCF, starting with the Ito stochastic equations for stochastic Grassmann
phase space variables. The stochastic averages that determine the matrices
that relate the stochastic averages for QCF at a later time (or temperature)
to the stochastic averages at an earlier time have been carried out
analytically, leading to standard non-stochastic matrices. The latter matrices
can then be used to study numerically the evolution of the QCF over a finite
time or temperature interval. Such calculations will be presented in a
subsequent paper.

The plan of the paper is as follows. In Section 3 we present the basic
equations for GPST in the context of the BEC/BCS crossover both for time an
temperature evolution, stating from the Hamiltonian, proceding via the FFPE
for the Grassmann distribution functional and ending with the Ito SFE both for
stochastic position fields and stochastic momentum fields. In Section 4 we
present the general expressions for the QCF for a single Cooper pair and for
two Cooper pairs. In Section 5 we present the results for the temperature and
time evolution for the QCF for a single Cooper pair, correct to first order
and with the stochastic averages carried out analytically. In Section 6, we
present the analogous results for the QCF describing two Cooper pairs.
Equations for carrying out the time or temperature evolution over a finite
interval based on the matrix $M$ obtained in the previous two Sections are set
out in Section 7. Finally, a Summary and Conclusions from the work is
presented in Section 8. Detailed derivations are set out in Appendices 10 and
11.

##  3 Grassman Phase Space Theory for Fields

In this Section we set out the basic equations for applying Grassmann phase
space theory to the topic of the BEC/BCS crossover in ultracold atomic fermi
gases. The general methodology and a full description of Grassmann algebra and
calculus is set out in the following textbook [19] and articles [27], [28],
[29] . We begin by writing the Hamiltonian involved in terms of fermion field
operators, then introduce the quantum density operator and its $B$
representation distribution functional involving fermion Grassmann field
functions, followed by the two basic equations the density operator satisfies
for describing time and temperature evolution. The quantum correlation
functions (QCF) relevant to fermion position measurements for the size of and
correlation between Cooper pairs are then defined. We set out the functional
Fokker-Planck equations (FFPE) for the distribution functional in the two
types of evolution, and then introduce the corresponding Ito stochastic field
equations that are equivalent to the FFPE. We finally obtain the Ito
stochastic equations for momentum fields in the case of untrapped fermions.

###  3.1 Hamiltonian and Field Operators

We consider a trapped Fermi gas of spin $1/2$ fermionic atoms with spin
conserving collisions of zero range between pairs of atoms .Here there are two
distinct internal states corresponding to spin up and spin down atoms,
designated by $\alpha=u(\uparrow),d(\downarrow)$, with $-\alpha$ referring to
the opposite spin state. The Fermi gas is isolated from the environment so no
relaxation effects are involved. This model was also considered by Plimak et
al [26].

The Hamiltonian is written in terms of the field operators
$\hat{\Psi}_{\alpha}(r)$, $\hat{\Psi}_{\alpha}(r)^{{\dagger}}$ as

$\displaystyle\hat{H}_{f}\mbox{\rule{-1.42262pt}{0.0pt}}$ $\displaystyle=$
$\displaystyle\mbox{\rule{-2.84526pt}{0.0pt}}\int\mbox{\rule{-2.84526pt}{0.0pt}}dr\left(\sum_{\alpha}\frac{\hbar^{2}}{2m}\nabla\hat{\Psi}_{\alpha}(r)^{{\dagger}}\mbox{\rule{-1.42262pt}{0.0pt}}\cdot\mbox{\rule{-1.42262pt}{0.0pt}}\nabla\hat{\Psi}_{\alpha}(\
r)+\sum_{\alpha}\hat{\Psi}_{\alpha}(r)^{{\dagger}}V_{\alpha}\hat{\Psi}_{\alpha}(r)\right.$
(1) $\displaystyle\hskip
56.9055pt\left.+\frac{g_{f}}{2}\sum_{\alpha}\hat{\Psi}_{\alpha}(r)^{{\dagger}}\hat{\Psi}_{-\alpha}(r)^{{\dagger}}\hat{\Psi}_{-\alpha}(r)\hat{\Psi}_{\alpha}(r)\right)$
$\displaystyle=$ $\displaystyle\widehat{{\small K}}+\widehat{{\small
V}}+\widehat{{\small U}}$ (2)

and is the sum of kinetic energy, trap potential energy and collision
interaction energy terms. Here $m$ is the mass of the fermionic atom,
$V_{\alpha}$ is the trap potential for spin $\alpha$ and $g$ is the fermion-
fermion interaction constant - which is the related to the $s$-wave scattering
amplitude $a_{s}$ via $g=4\pi a_{s}\hbar^{2}/m$. For simplicity $r$ refers to
position in 3D. The Hamiltonian commutes with the number operator
$\widehat{N}$, where

$\widehat{N}=\sum_{\alpha}\int\mbox{\rule{-2.84526pt}{0.0pt}}dr\,\left(\hat{\Psi}_{\alpha}(r)^{{\dagger}}\,\hat{\Psi}_{\alpha}(r)\right)$
(3)

###  3.2 B Distribution Functional

For fermionic systems with spin components $\alpha$ the field operators
$\hat{\Psi}_{\alpha}(r)$, $\hat{\Psi}_{\alpha}(r)^{{\dagger}}$ may be expanded
in a complete orthonormal set of single fermion mode functions $\phi_{\alpha
i}(r)$ for each component $\alpha$ as

$\widehat{\psi}_{\alpha}(r)=\mathop{\displaystyle\sum}\limits_{i}\widehat{c}_{\alpha
i}\;\phi_{\alpha
i}(r)\qquad\hat{\Psi}_{\alpha}(r)^{{\dagger}}=\mathop{\displaystyle\sum}\limits_{i}\widehat{c}_{\alpha
i}^{{\dagger}}\;\phi_{\alpha i}^{\ast}(r)$ (4)

where $\widehat{c}_{\alpha i}$, $\widehat{c}_{\alpha i}^{{\dagger}}$ are
fermionic annihilation, creation operators for the mode $\alpha i$. The
orthogonality and completeness properties are

$\int dr\,\phi_{\alpha i}^{\ast}(r)\phi_{\alpha
j}(r)=\delta_{ij}\qquad\sum_{i}\phi_{\alpha i}(r)\phi_{\alpha
i}^{\ast}(s)=\delta(r-s)$ (5)

In Grassmann phase space field theory a quantum density operator
$\widehat{\rho}$ can be represented by its $B$ distribution functional $B[\
\psi(r)]$ in which the field operators $\hat{\Psi}_{\alpha}(r)$,
$\hat{\Psi}_{\alpha}(r)^{{\dagger}}$ are represented via Grassmann fields
$\psi_{\alpha}(r)$ and $\psi_{\alpha}^{+}(r)$ respectively, where
$\psi(r)\equiv\\{\psi_{\alpha}(r),\psi_{\alpha}^{+}(r)\\}$. Here the Grassmann
fields are expanded as

$\psi_{\alpha}(r)=\mathop{\displaystyle\sum}\limits_{i}g_{\alpha
i}\;\phi_{\alpha
i}(r)\qquad\psi_{\alpha}^{+}(r)=\mathop{\displaystyle\sum}\limits_{i}g_{\alpha
i}^{+}\;\phi_{\alpha i}^{\ast}(r)$ (6)

where $g_{\alpha i}$, $g_{\alpha i}^{+}$ are Grassmann phase space variables
for the mode $\alpha i$.

The canonical form of a density operator for the $B$ distribution is written
in terms of Bargmann coherent state projectors $\widehat{\Lambda}[\psi(r)]$ in
terms of a Grassmann functional integral as

$\widehat{\rho}=\mathop{\displaystyle\int}\mathop{\displaystyle\prod}\limits_{\alpha}d\psi_{\alpha}^{+}\,d\psi_{\alpha}\,B[\
\psi(r)]\,\widehat{\Lambda}[\psi(r)]$ (7)

where

$\widehat{\Lambda}[\psi(r)]=\left|\psi(r)\right\rangle_{B}\;\left\langle\psi^{+}(r)^{\ast}\right|_{B}$
(8)

and the Bargmann states are given by

$\displaystyle\left|\psi(r)\right\rangle_{B}$ $\displaystyle=$
$\displaystyle\exp(\mathop{\displaystyle\prod}\limits_{\alpha}\mathop{\displaystyle\sum}\limits_{i}\widehat{c}_{\alpha
i}^{{\dagger}}\;g_{\alpha
i})\;\left|0\right\rangle=\mathop{\displaystyle\prod}\limits_{\alpha
i}(\left|0_{\alpha i}\right\rangle+g_{\alpha i}\left|1_{\alpha
i}\right\rangle)$ $\displaystyle\left|\psi^{+}(r)^{\ast}\right\rangle_{B}$
$\displaystyle=$
$\displaystyle\exp(\mathop{\displaystyle\prod}\limits_{\alpha}\mathop{\displaystyle\sum}\limits_{i}\widehat{c}_{\alpha
i}^{{\dagger}}\;(g_{\alpha
i}^{+})^{\ast})\;\left|0\right\rangle=\mathop{\displaystyle\prod}\limits_{\alpha
i}(\left|0_{\alpha i}\right\rangle+(g_{\alpha i}^{+})^{\ast}\left|1_{\alpha
i}\right\rangle)$ (9)

with $\left|0\right\rangle$ being the vacuum state with no fermions and
$\left|0_{\alpha i}\right\rangle$, $\left|1_{\alpha i}\right\rangle$ being one
fermion states with $0$, $1$ fermion in the $\alpha i$ mode.

The Bargmann states are eigenstates for the field operators

$\widehat{\psi}_{\alpha}(r)\,\left|\psi(r)\right\rangle_{B}=\psi_{\alpha}(r)\,\left|\psi(r)\right\rangle_{B}\qquad{}_{B}\left\langle\psi^{+}(r)^{\ast}\right|\,\hat{\Psi}_{\alpha}(r)^{{\dagger}}=_{B}\left\langle\psi^{+}(r)^{\ast}\right|\,\psi_{\alpha}^{+}(r)$
(10)

and are not normalised to unity. In fact

${}_{B}\left\langle\chi^{+}(r)^{\ast}|\,\psi(r)\right\rangle_{B}=\exp\mathop{\displaystyle\sum}\limits_{\alpha
i}h_{\alpha i}^{+}\,g_{\alpha i}$ (11)

where the Grassman field $\chi^{+}(r)$ involves phase space amplitudes
$h_{\alpha i}^{+}$. The Grassmann functional integral
$\mathop{\displaystyle\int}\mathop{\displaystyle\prod}\limits_{\alpha}d\psi_{\alpha}^{+}\,d\psi_{\alpha}\equiv\mathop{\displaystyle\int}d\psi^{+}\,d\psi$
involving the functional $B[\ \psi(r)]$ is equivalent to the multiple
Grassmann phase space integral over the phase space variables
$\mathop{\displaystyle\prod}\limits_{\alpha}dg_{\alpha 1}^{+}dg_{\alpha
2}^{+}....\,dg_{\alpha n}^{+}dg_{\alpha n}....dg_{\alpha 2}dg_{\alpha 1}$ of
the distribution function $B(g)$ (where $g\equiv\\{g_{\alpha i},g_{\alpha
i}^{+}\\}$) which is equivalent to the functional $B[\ \psi(r)]$.

###  3.3 Evolution of Density Operator

There are two cases of interest. The first is where the density operator
describes time evolution in a non-equilibrium system, in which case
$\widehat{\rho}$ satisfies the Liouville-von Neumann equation

$\frac{\partial}{\partial
t}\widehat{\rho}=\frac{-i}{\hbar}[\widehat{H}_{f}\,,\widehat{\rho}]$ (12)

involving the Hamiltonian $\widehat{H}_{f,}$ and the second where the density
operator describes temperature evolution in an equilibrium system, in which
case the density operator is given by the expression for the canonical
ensemble

$\displaystyle\widehat{\rho}$ $\displaystyle=$
$\displaystyle\exp(-\beta\widehat{H}_{f,})/Z$ (13) $\displaystyle Z$
$\displaystyle=$ $\displaystyle Tr\exp(-\beta\widehat{H}_{f})$ (14)

Here $Z$ is the partition function and $\beta=1/(k_{B}T)$, where $T$ is the
temperature. The numerator $\widehat{\sigma}=\exp(-\beta$ $\widehat{H}_{f,})$
satisfies a Matsubara equation.

$\frac{\partial}{\partial\beta}\widehat{\sigma}=\frac{-1}{2}[\widehat{H}_{f}\,,\widehat{\sigma}]_{+}$
(15)

In both cases $\widehat{\rho}$ and $\widehat{\sigma}$ can be represented via a
$B$ distribution, as in Eq.(7).

###  3.4 Quantum Correlation Functions and B Distribution Functional

Physical quantities of interest can be expressed as quantum correlation
functions involving the field operators $\hat{\Psi}_{\alpha}(r)$,
$\hat{\Psi}_{\alpha}(r)^{{\dagger}}$ and the density operator $\widehat{\rho}$
and $\widehat{\sigma}$. These QCF can then be written as Grassmann functional
integrals involving the $B$ distribution functional. Thus we define the QCF as

$\displaystyle
G^{(p,q)}(\alpha_{1}r_{1},\alpha_{2}r_{2}\cdots\alpha_{p}r_{p};\beta_{q}s_{q}\cdots\beta_{2}s_{2},\beta_{1}s_{1})$
(16) $\displaystyle=$
$\displaystyle\langle\hat{\Psi}_{\alpha_{1}}(r_{1})^{{\dagger}}\cdots\hat{\Psi}_{\alpha_{p}}(r_{p})^{{\dagger}}\hat{\Psi}_{\beta_{q}}(s_{q})\cdots\hat{\Psi}_{\beta_{1}}(s_{1})\rangle$
$\displaystyle=$
$\displaystyle\mbox{Tr}(\hat{\Psi}_{\beta_{q}}(s_{q})\cdots\hat{\Psi}_{\beta_{1}}(s_{1})\,\hat{\rho}\,\hat{\Psi}_{\alpha_{1}}(r_{1})^{{\dagger}}\cdots\hat{\Psi}_{\alpha_{p}}(r_{p})^{{\dagger}}).\quad$

and in terms of Grassmann functional integrals the QCF are

$\displaystyle
G^{(p,q)}(\alpha_{1}r_{1},\alpha_{2}r_{2}\cdots\alpha_{p}r_{p};\beta_{q}s_{q}\cdots\beta_{2}s_{2},\beta_{1}s_{1})$
$\displaystyle=$ $\displaystyle\int\int
d\psi^{+}d\psi\,\psi_{\beta_{q}}(s_{q})\cdots\psi_{\beta_{1}}(\
s_{1})B[\psi(r),\psi^{+}(r)]\psi_{\alpha_{1}}^{+}(r_{1})\cdots\psi_{\alpha_{p}}^{+}(r_{p}))$

for both the $\widehat{\rho}$ and $\widehat{\sigma}$ cases.

###  3.5 Functional Fokker-Planck Equations

The functional Fokker-Planck equation for the $B$ distribution functional $B[\
\psi(r)]$ is derived via using the correspondence rules in conjunction with
functional differentiation rules. We have for both the $\widehat{\rho}$ and
$\widehat{\sigma}$ cases

$\displaystyle\hat{\rho}$ $\displaystyle\Rightarrow$
$\displaystyle\widehat{{\small\Psi}}_{\alpha}(r)\,\hat{\rho}\quad
B[\psi]\Rightarrow\psi_{\alpha}(r)\,B[\psi]\qquad\qquad\hat{\rho}\Rightarrow\hat{\rho}\,\widehat{{\small\Psi}}_{\alpha}(r)\qquad
B[\psi]\Rightarrow
B[\psi]\,(+\frac{\overleftarrow{\delta}}{\delta\psi_{\alpha}^{+}(r)})$
$\displaystyle\hat{\rho}$ $\displaystyle\Rightarrow$
$\displaystyle\widehat{{\small\Psi}}_{\alpha}^{{\dagger}}(r)\,\hat{\rho}\quad
B[\psi]\Rightarrow(+\frac{\overrightarrow{\delta}}{\delta\psi_{\alpha}(r)})\,B[\psi]\qquad\hat{\rho}\Rightarrow\hat{\rho}\,\widehat{{\small\Psi}}_{\alpha}^{{\dagger}}(r)\quad
B[\psi]\Rightarrow\,B[\psi]\,\psi_{\alpha}^{+}(r)$ (18)

####  3.5.1 Time Evolution Case

In the case of the density operator $\widehat{\rho}$ satisfying the Liouville-
von Neumann equation the time evolution of the $B$ distribution functional is
as follows. There are terms corresponding to the kinetic, trap potential and
fermion-fermion interaction terms in the Hamiltonian. The latter term has been
written as a double space integral via the introduction of a Dirac delta
function to comply with the requirements for deriving the Ito stochastic
equations [29] from the general theory [19], [27], [28] .

The kinetic energy term is

$\displaystyle\left(\frac{\partial}{\partial t}B[\psi(r)]\right)_{K}$
$\displaystyle=$ $\displaystyle\frac{i}{\hbar}\int
ds\,\left[\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{u}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{d}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}(s)}\right\\}\right.$
$\displaystyle\left.-\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{u}^{+}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(s)}\right\\}\right.$
$\displaystyle\left.-\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{d}^{+}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(s)}\right\\}\right]$

and only contributes to the drift term.

The trap potential energy term is

$\displaystyle\left(\frac{\partial}{\partial t}B[\psi(r)]\right)_{V}$
$\displaystyle=$ $\displaystyle\frac{i}{\hbar}\int{\small
ds}\,\left[-\left\\{\left(V_{u}\psi_{u}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}(s)}\right\\}-\left\\{\left(V_{d}\psi_{d}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\left(V_{u}\psi_{u}^{+}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(s)}\right\\}+\left\\{\left(V_{d}\psi_{d}^{+}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(s)}\right\\}\right]$

and also only contributes to the drift term.

The fermion-fermion interaction term is

$\displaystyle\left(\frac{\partial}{\partial t}B[\psi(r)]\right)_{U}$ (21)
$\displaystyle{\small=}$
$\displaystyle\frac{i}{\hbar}\frac{g}{2}\int\int{\small
ds}\,dr\,\left[\left\\{\psi_{d}(s)\,\delta(r-s)\,\psi_{u}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{d}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{u}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\psi_{u}(s)\,\delta(r-s)\,\psi_{d}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{u}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{d}(s)}\right\\}\right.$
$\displaystyle\left.-\left\\{\psi_{d}^{+}(s)\,\delta(r-s)\,\psi_{u}^{+}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(s)}\right\\}\right.$
$\displaystyle\left.-\left\\{\psi_{u}^{+}(s)\,\delta(r-s)\,\psi_{d}^{+}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(s)}\right\\}\right]$

and is the only contribution to the diffusion term.

####  3.5.2 Temperature Evolution Case

In the case of the density operator $\widehat{\sigma}$ satisfying the
Matsubara equation the temperature evolution of the $B$ distribution
functional is as follows.

The kinetic energy term is

$\displaystyle\left(\frac{\partial}{\partial\beta}B[\psi(r)]\right)_{K}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\int
ds\,\left[\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{u}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{d}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{u}^{+}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\left(\frac{\hbar^{2}}{2m}\nabla^{2}{\small\psi}_{d}^{+}{\small(s)\,}B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(s)}\right\\}\right]$

and only contributes to the drift term.

The trap potential energy term is

$\displaystyle\left(\frac{\partial}{\partial\beta}B[\psi(r)]\right)_{V}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\int{\small
ds}\,\left[+\left\\{\left(V_{u}\psi_{u}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}(s)}\right\\}+\left\\{\left(V_{d}\psi_{d}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\left(V_{u}\psi_{u}^{+}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(s)}\right\\}+\left\\{\left(V_{d}\psi_{d}^{+}(s)\,B[\
\psi(r)]\right)\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(s)}\right\\}\right]$

and also only contributes to the drift term.

The fermion-fermion interaction term is

$\displaystyle\left(\frac{\partial}{\partial\beta}B[\psi(r)]\right)_{U}$ (24)
$\displaystyle{\small=}$ $\displaystyle{\small+}\frac{g}{4}\int\int{\small
ds}\,dr\,\left[\left\\{\psi_{d}(s)\,\delta(r-s)\,\psi_{u}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{d}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{u}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\psi_{u}(s)\,\delta(r-s)\,\psi_{d}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{u}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{d}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\psi_{d}^{+}(s)\,\delta(r-s)\,\psi_{u}^{+}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(s)}\right\\}\right.$
$\displaystyle\left.+\left\\{\psi_{u}^{+}(s)\,\delta(r-s)\,\psi_{d}^{+}(r)\,B[\psi(r)]\frac{\overleftarrow{\delta}}{\delta\psi_{u}^{+}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{d}^{+}(s)}\right\\}\right]$

and is the only contribution to the diffusion term. Note the differences
between the Liouville-von Neumann and Matusubara cases. Apart from the overall
multiplier being related to $-i/\hbar$ rather than $-1/2$, the pairs of terms
involving $\psi_{\alpha}(r)$ and the pairs of terms involving
$\psi_{\alpha}^{+}(r)$ have the same sign in the Matsubara case and opposite
signs in the other case.

###  3.6 Ito Stochastic Position Field Equations

The Ito stochastic field equations are derived from the FFPE via using the
first order functional derivative (or drift) terms to give the classical (or
non-noise) terms in the Ito SFE and via using the second order (or diffusion)
terms to determine the noise (or fluctuation) terms in the Ito SFE. The latter
terms involve Gaussian-Markoff stochastic noise quantities and require the
Takagi factorisation of the diffusion matrix.

The basic principle is that the time (or temperature) evolution of the $B$
distribution functional is equivalent to the time (or temperature) evolution
of stochastic Grassmann fields
$\widetilde{\psi}_{\alpha}(r,t),\widetilde{\psi}_{\alpha}^{+}(r,t)$ (or
$\widetilde{\psi}_{\alpha}(r,\beta),\widetilde{\psi}_{\alpha}^{+}(r,\beta)$)
which now replace the Grassman phase space fields
$\psi_{\alpha}(r),\psi_{\alpha}^{+}(r)$, and where the evolution is now
described by Ito stochastic field equations. These are chosen so that the
stochastic averages of the
$\widetilde{\psi}_{\alpha}(r,t),\widetilde{\psi}_{\alpha}^{+}(r,t)$ (or
$\widetilde{\psi}_{\alpha}(r,\beta),\widetilde{\psi}_{\alpha}^{+}(r,\beta)$)
gives the same result as the Grassmann functional integrals involving the
$\psi_{\alpha}(r),\psi_{\alpha}^{+}(r)$. In particular, for the time evolution
case the QCF are now given by

$\displaystyle
G^{(p,q)}(\alpha_{1}r_{1},\alpha_{2}r_{2}\cdots\alpha_{p}r_{p};\beta_{q}s_{q}\cdots\beta_{2}s_{2},\beta_{1}s_{1})$
(25) $\displaystyle=$
$\displaystyle\mbox{Tr}(\hat{\Psi}_{\beta_{q}}(s_{q})\cdots\hat{\Psi}_{\beta_{1}}(s_{1})\,\hat{\rho}\,\hat{\Psi}_{\alpha_{1}}(r_{1})^{{\dagger}}\cdots\hat{\Psi}_{\alpha_{p}}(r_{p})^{{\dagger}}).\quad$
$\displaystyle=$ $\displaystyle\int\int
d\psi^{+}d\psi\,\psi_{\beta_{q}}(s_{q})\cdots\psi_{\beta_{1}}(\
s_{1})B[\psi(r),\psi^{+}(r)]\psi_{\alpha_{1}}^{+}(r_{1})\cdots\psi_{\alpha_{p}}^{+}(r_{p}))$
$\displaystyle=$
$\displaystyle\overline{\widetilde{\psi}_{\beta_{q}}(s_{q},t)\cdots\widetilde{\psi}_{\beta_{1}}(\
s_{1},t)\;\widetilde{\psi}_{\alpha_{1}}^{+}(r_{1},t)\cdots\widetilde{\psi}_{\alpha_{p}}^{+}(r_{p},t)}$

where the bar indicates a stochastic average. The same expressions apply for
the temperature evolution case ,but with $t$ being replaced by $\beta$.

To set out the FFPE and the related Ito SFE it is useful to modify the
notation for the field functions. The indices $A,B$ are inroduced where each
takes on the values $1,2$ and where (which $\alpha,\beta$ specifying $u$ or
$d$ as before) $\psi_{\alpha 1}(s)\equiv\psi_{\alpha}(s)$ and $\psi_{\alpha
2}(s)\equiv\psi_{\alpha}^{+}(s)$. Similarly, $\widetilde{\psi}_{\alpha
1}(s,t)\equiv\widetilde{\psi}_{\alpha}(s,t)$ and $\widetilde{\psi}_{\alpha
2}(s,t)\equiv\widetilde{\psi}_{\alpha}^{+}(s,t)$ for the stochastic fields.

For the case of time evolution we may write the functional Fokker-Planck
equation in the general form

$\displaystyle\frac{\partial}{\partial t}B[\psi]$ $\displaystyle=$
$\displaystyle-\mathop{\textstyle\sum}\limits_{\alpha
A}\mathop{\displaystyle\int}dr\,(A_{\alpha
A}[\psi(r),r]\,B[\psi])\frac{\overleftarrow{\delta}}{\delta\psi_{\alpha
A}(r)}$
$\displaystyle{\small+}\frac{1}{2}\mathop{\textstyle\sum}\limits_{\alpha
A,\beta B}\mathop{\displaystyle\iint}ds\,dr\,(D_{\alpha A\,\beta
B}[\psi(s),s;\psi(r),r]\,B[\psi])\frac{\overleftarrow{\delta}}{\delta\psi_{\beta
B}(r)}\frac{\overleftarrow{\delta}}{\delta\psi_{\alpha A}(s)}$

and then the Ito stochastic field equations can be written as the sum of a
classical and a noise contribution for the change in $\widetilde{\psi}_{\alpha
A}$ as

$\displaystyle\delta\widetilde{\psi}_{\alpha A}(r,t)$ $\displaystyle\equiv$
$\displaystyle\widetilde{\psi}_{\alpha A}(r,t+\delta
t)-\widetilde{\psi}_{\alpha A}(r,t)$ (27) $\displaystyle=$ $\displaystyle-
A_{\alpha A}[\widetilde{\psi}(r,t),r]\,\delta
t+\mathop{\textstyle\sum}\limits_{a}B_{a}^{\alpha
A}[\widetilde{\psi}(r,t),r]\;\delta\omega_{a}(t_{+})$ $\displaystyle=$
$\displaystyle\left(\delta\widetilde{\psi}_{\alpha A}({\small
r})\right)_{class}+\left(\delta\widetilde{\psi}_{\alpha A}({\small
r})\right)_{noise}$ (28)

where the matrix $B$ is linked to the diffusion matrix $D$ via Takagi
factorisation as

$\displaystyle(BB^{T})_{\alpha A\,\beta B}$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{a}B_{a}^{\alpha
A}[\widetilde{\psi}(s,t),s]\,B_{a}^{\beta B}[\widetilde{\psi}(r,t),r]$ (29)
$\displaystyle=$ $\displaystyle D_{\alpha A\,\beta
B}[\widetilde{\psi}(s,t),s;\widetilde{\psi}(r,t),r]$

The drift and diffusion matrices can be read off from Eqs.
(LABEL:Eq.FFPEKineticFermiFldModel), (LABEL:Eq.FFPEPotentialFermiFldModel) and
(21). In determining the matrix $B$ the delta function
$\delta(r-s)=\frac{1}{V}\sum\limits_{q}\exp iq\cdot(r-s)$ expressed in terms
of box normalised plane waves, with $V=L^{3}$ and (for 3D)
$q\equiv\\{q_{x}(\frac{2\pi}{L}),q_{y}(\frac{2\pi}{L}),q_{z}(\frac{2\pi}{L})\\}$
\- where $q_{x},q_{y},q_{z}$ are all integers.is used in the diffusion matrix.
The details for determining the matrix $B$ are set out in Appendix 10.

The Wiener increments $\delta\omega_{a}(t_{+})$ are defined in terms of
Gausian-Markoff random noise terms $\Gamma_{a}(t)$ as

$\delta\omega_{a}(t_{+})=\int_{t}^{t+\delta t}dt_{1}\,\Gamma_{a}(t_{1})$ (30)

Properties of the Gaussian-Markoff random noise terms are set out in Appendix
10. From these results we can then show that the Wiener increments have the
following properties

$\displaystyle\overline{\delta\omega_{{\small a}}(t_{+})}\,$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\overline{\delta\omega_{{\small
a}}(t_{+})\delta\omega_{b}(t_{+})}\,$ $\displaystyle=$
$\displaystyle\delta_{{\small ab}}\;\delta t$
$\displaystyle\overline{\delta\omega_{a}(t_{+})\delta\omega_{b}(t_{+})\delta\omega_{c}(t_{+})}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\overline{\delta\omega_{a}(t_{+})\delta\omega_{b}(t_{+})\delta\omega_{c}(t_{+})\delta\omega_{d}(t_{+})}$
$\displaystyle=$ $\displaystyle\overline{\delta\omega_{{\small
a}}(t_{+})\delta\omega_{b}(t_{+})}\;\overline{\delta\omega_{{\small
c}}(t_{+})\delta\omega_{d}(t_{+})}+\overline{\delta\omega_{{\small
a}}(t_{+})\delta\omega_{c}(t_{+})}\;\,\overline{\delta\omega_{{\small
b}}(t_{+})\delta\omega_{d}(t_{+})}$
$\displaystyle+\overline{\delta\omega_{{\small
a}}(t_{+})\delta\omega_{d}(t_{+})}\;\,\overline{\delta\omega_{{\small
b}}(t_{+})\delta\omega_{c}(t_{+})}$

showing that the stochastic averages of a single $\delta\omega_{a}$ is zero
and the stochastic average of the product of two $\delta\omega_{a}$’s is zero
if they are different and equal to $\delta t$ if they are the same. In
addition, the stochastic averages of products of odd numbers of
$\delta\omega_{a}$ are zero and stochastic averages of products of even
numbers of $\delta\omega_{a}$ are the sums of products of stochastic averages
of pairs of $\delta\omega_{a}$.

For the case of temperature evolution, the same forms apply but with $t$ being
replaced by $\beta$.

####  3.6.1 Time Evolution Case

For the time evolution case the stochastic fields at time $t+\delta t$ can be
related to the stochastic fields at time $t$ as follows. We have now reverted
to the original notation.

We have (see Sect. 10.2 for details)

$\displaystyle\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\psi}_{u}(s,t+\delta
t)$\\\ \hline\cr$\widetilde{\psi}_{d}(s,t+\delta t)$\\\
\hline\cr$\widetilde{\psi}_{u}^{+}(s,t+\delta t)$\\\
\hline\cr$\widetilde{\psi}_{d}^{+}(s,t+\delta t)$\\\
\hline\cr\end{tabular}\right]$ (36) $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|l|l}\hline\cr$\Theta_{u,u}(\delta
t)$&$\Theta_{u,d}(\delta t)$&$0$&$0$\\\ \hline\cr$\Theta_{d,u}(\delta
t)$&$\Theta_{d,d}(\delta t)$&$0$&$0$\\\
\hline\cr$0$&$0$&$\Theta_{u,u}^{+}(\delta t)$&$\Theta_{u,d}^{+}(\delta t)$\\\
\hline\cr$0$&$0$&$\Theta_{d,u}^{+}(\delta t)$&$\Theta_{d,d}^{+}(\delta t)$\\\
\hline\cr\end{tabular}\right]\;\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\psi}_{u}(s,t)$\\\
\hline\cr$\widetilde{\psi}_{d}(s,t)$\\\
\hline\cr$\widetilde{\psi}_{u}^{+}(s,t)$\\\
\hline\cr$\widetilde{\psi}_{d}^{+}(s,t)$\\\ \hline\cr\end{tabular}\right]$
(45)

where the $\Theta$ quantities are

$\displaystyle\left[\Theta(\delta t)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$1-\frac{i}{\hbar}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{u}\right\\}\delta
t$&$+\lambda\sum\limits_{q}\exp(+iqs)\left\\{\delta\widetilde{\omega}_{u,d}^{q}+i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$\\\
\hline\cr$+\lambda\sum\limits_{q}\exp(-iqs)\left\\{\delta\widetilde{\omega}_{u,d}^{q}-i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$&$1-\frac{i}{\hbar}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{d}\right\\}\delta
t$\\\ \hline\cr\end{tabular}\right]$ (49)
$\displaystyle\left[\Theta^{+}(\delta t)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$1+\frac{i}{\hbar}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{u}\right\\}\delta
t$&$+\lambda\sum\limits_{q}\exp(-iqs)\left\\{\delta\widetilde{\omega}_{u+,d+}^{q}+i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$\\\
\hline\cr$+\lambda\sum\limits_{q}\exp(+iqs)\left\\{-\delta\widetilde{\omega}_{u+,d+}^{q}+i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$&$1+\frac{i}{\hbar}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{d}\right\\}\delta
t$\\\ \hline\cr\end{tabular}\right]$ (52)

where $\lambda=\sqrt{\frac{ig}{2\hslash V}}$. These quantities involve the
Laplacian, the trap potentials, the time interval $\delta t$ as well as four
stochastic Wiener increments (30) $\delta\widetilde{\omega}_{u,d}^{q}$,
$\delta\widetilde{\omega}_{d,u}^{q}$, $\delta\widetilde{\omega}_{u+,d+}^{q}$
and $\widetilde{\omega}_{d+,u+}^{q}$ for each $q$. The latter satisfy
conditions analogous to (LABEL:Eq.PropsWienerIncrem). Note the sum over $q$
and that there is no change in the position coordinate $s$. This is a
consequence of the zero range interaction.

####  3.6.2 Temperature Evolution Case

For the temperature evolution case the stochastic fields at temperature
$\beta+\delta\beta$ can be related to the stochastic fields at temperature
$\beta$ as follows. We have now reverted to the original notation.

We have (see Sect. 10.3 for details)

$\displaystyle\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\psi}_{u}(s,\beta+\delta\beta)$\\\
\hline\cr$\widetilde{\psi}_{d}(s,\beta+\delta\beta)$\\\
\hline\cr$\widetilde{\psi}_{u}^{+}(s,\beta+\delta\beta)$\\\
\hline\cr$\widetilde{\psi}_{d}^{+}(s,\beta+\delta\beta)$\\\
\hline\cr\end{tabular}\right]$ (58) $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|l|l}\hline\cr$\Theta_{u,u}(\delta\beta)$&$\Theta_{u,d}(\delta\beta)$&$0$&$0$\\\
\hline\cr$\Theta_{d,u}(\delta\beta)$&$\Theta_{d,d}(\delta\beta)$&$0$&$0$\\\
\hline\cr$0$&$0$&$\Theta_{u,u}^{+}(\delta\beta)$&$\Theta_{u,d}^{+}(\delta\beta)$\\\
\hline\cr$0$&$0$&$\Theta_{d,u}^{+}(\delta\beta)$&$\Theta_{d,d}^{+}(\delta\beta)$\\\
\hline\cr\end{tabular}\right]\;\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\psi}_{u}(s,\beta)$\\\
\hline\cr$\widetilde{\psi}_{d}(s,\beta)$\\\
\hline\cr$\widetilde{\psi}_{u}^{+}(s,\beta)$\\\
\hline\cr$\widetilde{\psi}_{d}^{+}(s,\beta)$\\\ \hline\cr\end{tabular}\right]$
(67)

where the $\Theta$ quantities are

$\displaystyle\left[\Theta(\delta\beta)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$1-\frac{1}{2}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{u}\right\\}\delta\beta$&$+\eta\sum\limits_{q}\exp(+iqs)\left\\{\delta\widetilde{\omega}_{u,d}^{q}+i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$\\\
\hline\cr$+\eta\sum\limits_{q}\exp(-iqs)\left\\{\delta\widetilde{\omega}_{u,d}^{q}-i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$&$1-\frac{1}{2}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{d}\right\\}\delta\beta$\\\
\hline\cr\end{tabular}\right]$ (71)
$\displaystyle\left[\Theta^{+}(\delta\beta)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$1-\frac{1}{2}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{u}\right\\}\delta\beta$&$+\eta\sum\limits_{q}\exp(-iqs)\left\\{\delta\widetilde{\omega}_{u+,d+}^{q}+i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$\\\
\hline\cr$+\eta\sum\limits_{q}\exp(+iqs)\left\\{\delta\widetilde{\omega}_{u+,d+}^{q}-i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$&$1-\frac{1}{2}\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{d}\right\\}\delta\beta$\\\
\hline\cr\end{tabular}\right]$ (74)

where $\eta=\sqrt{\frac{g}{8V}}$. These quantities involve the Laplacian, the
trap potentials, the temperature interval $\delta\beta$ as well as four
stochastic Wiener increments (30) $\delta\widetilde{\omega}_{u,d}^{q}$,
$\delta\widetilde{\omega}_{d,u}^{q}$, $\delta\widetilde{\omega}_{u+,d+}^{q}$
and $\widetilde{\omega}_{d+,u+}^{q}$ for each $q$. The latter satisfy
conditions analogous to (LABEL:Eq.PropsWienerIncrem). Note the sum over $q$
and that there is no change in the position coordinate $s$. This is a
consequence of the zero range interaction. Note the sign differences for both
the classical and noise terms between the time and temperature evolution
cases. There is also a change in the constants for these terms.

###  3.7 Ito Stochastic Momentum Field Equations

For the free field situation where $V_{u}=V_{d}=0$, it is convenient to
introduce stochastic momentum fields
$\tilde{\phi}_{\alpha}(k),\,\tilde{\phi}_{\alpha}^{+}(k)$, which are related
to the stochastic position fields
$\tilde{\psi}_{\alpha}(s),\tilde{\psi}_{\alpha}^{+}(s)$ via spatial Fourier
transforms involving box normalised plane waves $\exp(ik\cdot s)/\sqrt{V}$.
The wave numbers $k$ are discrete and given by
$k\equiv\\{k_{x}(\frac{2\pi}{L}),k_{y}(\frac{2\pi}{L}),k_{z}(\frac{2\pi}{L})\\}$
\- where $k_{x},k_{y},k_{z}$ are all integers and $V=L^{3}$. For ease of
notation, we leave the $t$ or $\beta$ dependence implicit.

$\displaystyle\tilde{\psi}_{\alpha}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{V}}\mathop{\displaystyle\sum}\limits_{k}\,\exp(ik\cdot
s)\,\tilde{\phi}_{\alpha}(k)$ $\displaystyle\tilde{\psi}_{\alpha}^{+}(s)$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{V}}\mathop{\displaystyle\sum}\limits_{k}\,\exp(-ik\cdot
s)\,\tilde{\phi}_{\alpha}^{+}(k)$ (76)

The inverse transformation is given as Eq. (198) in Appendix 10. The use of
stochastic momentum fields enables the $\nabla^{2}$ terms to be replaced by
$k^{2}$ when the stochastic position fields are replaced by stochastic
momentum fields in Eqs. (LABEL:Eq.ItoSFETimeEvoln) and
(LABEL:Eq.ItoSFETempEvn).

####  3.7.1 Time Evolution Case - Momentum Fields

For the time evolution case the stochastic momentum fields at time $t+\delta
t$ can be related to the stochastic momentum fields at time $t$ as follows. We
have now reverted to the original notation. To obtain this result we
substitute the expressions (76) for the stochastic position fields into Eq.
(LABEL:Eq.ItoSFETimeEvoln).

We have

$\displaystyle\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\phi}_{u}(k,t+\delta
t)$\\\ \hline\cr$\widetilde{\phi}_{d}(k,t+\delta t)$\\\
\hline\cr$\widetilde{\phi}_{u}^{+}(k,t+\delta t)$\\\
\hline\cr$\widetilde{\phi}_{d}^{+}(k,t+\delta t)$\\\
\hline\cr\end{tabular}\right]$ (81) $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|l|l}\hline\cr$F_{u,u}(k,l;q,\delta
t)$&$F_{u,d}(k,l;q,\delta t)$&$0$&$0$\\\ \hline\cr$F_{d,u}(k,l;q,\delta
t)$&$F_{d,d}(k,l;q,\delta t)$&$0$&$0$\\\
\hline\cr$0$&$0$&$F_{u,u}^{+}(k,l;q,\delta t)$&$F_{u,d}^{+}(k,l;q,\delta
t)$\\\ \hline\cr$0$&$0$&$F_{d,u}^{+}(k,l;q,\delta
t)$&$F_{d,d}^{+}(k,l;q,\delta t)$\\\
\hline\cr\end{tabular}\right]\;\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\phi}_{u}(l,t)$\\\
\hline\cr$\widetilde{\phi}_{d}(l,t)$\\\
\hline\cr$\widetilde{\phi}_{u}^{+}(l,t)$\\\
\hline\cr$\widetilde{\phi}_{d}^{+}(l,t)$\\\ \hline\cr\end{tabular}\right]$
(90)

where the $F$ quantities are

$\displaystyle\left[F(k,l;q,\delta t)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$\delta_{q,0}\;\delta_{k,l}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta
t\right)$&$+\lambda\;\delta_{(k-q),l}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q}+i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$\\\
\hline\cr$+\lambda\;\delta_{(k+q),l}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q}-i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$&$\delta_{q,0}\;\delta_{k,l}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta
t\right)$\\\ \hline\cr\end{tabular}\right]$ (94)
$\displaystyle\left[F^{+}(k,l;q,\delta t)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$\delta_{q,0}\;\delta_{k,l}\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta
t\right)$&$+\lambda\;\delta_{(k+q),l}\;\left\\{\delta\widetilde{\omega}_{u+,d+}^{q}+i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$\\\
\hline\cr$+\lambda\;\delta_{(k-q),l}\;\left\\{-\delta\widetilde{\omega}_{u+,d+}^{q}+i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$&$\delta_{q,0}\;\delta_{k,l}\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta
t\right)$\\\ \hline\cr\end{tabular}\right]$ (97)

where $\lambda=\sqrt{\frac{ig}{2\hbar V}}$. In each row of the right side of
Eq.(LABEL:Eq.ItoSFETimeEvoln) the sum over repeated indices $q,l$ is assumed.

####  3.7.2 Temperature Evolution Case - Momentum Fields

For the temperature evolution case the stochastic momentum fields at
temperature $\beta+\delta\beta$ can be related to the stochastic momentum
fields at temperature $\beta$ as follows. We have now reverted to the original
notation. To obtain this result we substitute the expressions (76) for the
stochastic position fields into Eq. (LABEL:Eq.ItoSFETempEvn).

We have

$\displaystyle\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\phi}_{u}(k,\beta+\delta\beta)$\\\
\hline\cr$\widetilde{\phi}_{d}(k,\beta+\delta\beta)$\\\
\hline\cr$\widetilde{\phi}_{u}^{+}(k,\beta+\delta\beta)$\\\
\hline\cr$\widetilde{\phi}_{d}^{+}(k,\beta+\delta\beta)$\\\
\hline\cr\end{tabular}\right]$ (103) $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|l|l}\hline\cr$F_{u,u}(k,l;q,\delta\beta)$&$F_{u,d}(k,l;q,\delta\beta)$&$0$&$0$\\\
\hline\cr$F_{d,u}(k,l;q,\delta\beta)$&$F_{d,d}(k,l;q,\delta\beta)$&$0$&$0$\\\
\hline\cr$0$&$0$&$F_{u,u}^{+}(k,l;q,\delta\beta)$&$F_{u,d}^{+}(k,l;q,\delta\beta)$\\\
\hline\cr$0$&$0$&$F_{d,u}^{+}(k,l;q,\delta\beta)$&$F_{d,d}^{+}(k,l;q,\delta\beta)$\\\
\hline\cr\end{tabular}\right]\;\left[\begin{tabular}[]{l}\hline\cr$\widetilde{\phi}_{u}(l,\beta)$\\\
\hline\cr$\widetilde{\phi}_{d}(l,\beta)$\\\
\hline\cr$\widetilde{\phi}_{u}^{+}(l,\beta)$\\\
\hline\cr$\widetilde{\phi}_{d}^{+}(l,\beta)$\\\ \hline\cr\end{tabular}\right]$
(112)

where the $F$ quantities are

$\displaystyle\left[F(k,l;q,\delta\beta)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$\delta_{q,0}\;\delta_{k,l}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta\beta\right)$&$+\eta\;\delta_{(k-q),l}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q}+i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$\\\
\hline\cr$+\eta\;\delta_{(k+q),l}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q}-i\delta\widetilde{\omega}_{d,u}^{q}\right\\}$&$\delta_{q,0}\;\delta_{k,l}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta\beta\right)$\\\
\hline\cr\end{tabular}\right]$ (116)
$\displaystyle\left[F^{+}(k,l;q,\delta\beta)\right]$ $\displaystyle=$
$\displaystyle\left[\begin{tabular}[]{|l|l|}\hline\cr$\delta_{q,0}\;\delta_{k,l}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta\beta\right)$&$+\eta\;\delta_{(k+q),l}\;\left\\{\delta\widetilde{\omega}_{u+,d+}^{q}+i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$\\\
\hline\cr$+\eta\;\delta_{(k-q),l}\;\left\\{\delta\widetilde{\omega}_{u+,d+}^{q}-i\delta\widetilde{\omega}_{d+,u+}^{q}\right\\}$&$\delta_{q,0}\;\delta_{k,l}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k^{2}}{2m}\right\\}\delta\beta\right)$\\\
\hline\cr\end{tabular}\right]$ (119)

where $\eta=\sqrt{\frac{g}{8V}}$. In each row of the right side of
Eq.(LABEL:Eq.ItoSMtmFTemper) the sum over repeated indices $q,l$ is assumed.

##  4 Quantum Correlation Functions for Cooper Pairs

In this Section we derive the expressions for the QCF’s involving in
determining the size of a single Cooper pair and determining the correlation
between the positions of the fermions in two Cooper pairs.

For the case of a single Cooper pair where the spin down fermion is at
position $r_{1}$ and the spin up fermion is at position $r_{2}$ the QCF that
can be used to describe the size of the Cooper pair is

$X(d\,r_{1},u\,r_{2})=\mbox{Tr}(\hat{\Psi}_{d}(r_{1})\,\hat{\Psi}_{u}(r_{2})\,\hat{\rho}\,\hat{\Psi}_{u}(r_{2})^{{\dagger}}\,\hat{\Psi}_{d}(r_{1})^{{\dagger}}).$
(121)

For the case of two Cooper pairs where the spin down fermions are at positions
$r_{1}$, $r_{2}$ and the spin up fermions are at positions $r_{3}$, $r_{4}$
the QCF that can be used to describe the correlation of the positions of the
fermions in the two Cooper pairs is

$\displaystyle X(d\,r_{1},d\,r_{2,}\,u\,r_{3},u\,r_{4})=$
$\displaystyle\mbox{Tr}(\hat{\Psi}_{d}(r_{1})\,\hat{\Psi}_{d}(r_{2})\,\hat{\Psi}_{u}(r_{3})\,\hat{\Psi}_{u}(r_{4})\,\hat{\rho}\,\hat{\Psi}_{u}(r_{4})^{{\dagger}}\,\hat{\Psi}_{u}(r_{3})^{{\dagger}}\,\hat{\Psi}_{d}(r_{2})^{{\dagger}}\,\hat{\Psi}_{d}(r_{1})^{{\dagger}}).$
(122)

###  4.1 One Cooper Pair: QCF

For the case of time evolution the stochastic expressions for
$X(d\,r_{1},u\,r_{2})$ can be written in terms of stochastic position or
momentum fields as

$\displaystyle X(d\,r_{1},u\,r_{2};t)$ $\displaystyle=$
$\displaystyle\overline{\widetilde{\psi}_{d}(r_{1},t)\,\widetilde{\psi}_{u}(r_{2},t)\;\widetilde{\psi}_{u}^{+}(r_{2},t)\,\widetilde{\psi}_{d}^{+}(r_{1},t)}$
(124) $\displaystyle=$
$\displaystyle\frac{1}{(\sqrt{V})^{4}}\mathop{\displaystyle\sum}\limits_{k_{1},k_{2},k_{3},k_{4}}\exp
i\\{(k_{1}-k_{4})\cdot r_{1}\\}\;\exp i\\{(k_{2}-k_{3})\cdot r_{2}\\}$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(k_{1},t)\,\widetilde{\phi}_{u}(k_{2},t)\;\widetilde{\phi}_{u}^{+}(k_{3},t)\,\widetilde{\phi}_{d}^{+}(k_{4},t)}$

The case of temperature evolution leads to the same forms with $t$ replaced by
$\beta$ in Eq.(124) and in Eq. (124), along with the additional factor $1/Z$
on the right side of Eq. (124) to allow for $\widehat{\rho}$ being replaced by
$\widehat{\sigma}$ when temperature evolution is involved.

The momentum stochastic correlation function of interest is

$X(d\,k_{1},u\,k_{2},u^{+}\,k_{3},d^{+}\,k_{4})=\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
(125)

where we have left the $t$ or $\beta$ dependences implicit to shorten the
notation.

###  4.2 Two Cooper Pairs: QCF

For the case of time evolution the stochastic expressions for
$X(d\,r_{1},d\,r_{2,}\,u\,r_{3},u\,r_{4})$ can be written in terms of
stochastic position or momentum fields as

$\displaystyle X(d\,r_{1},d\,r_{2,}\,u\,r_{3},u\,r_{4};t)$ $\displaystyle=$
$\displaystyle\overline{\widetilde{\psi}_{d}(r_{1},t)\,\widetilde{\psi}_{d}(r_{2},t)\,\widetilde{\psi}_{u}(r_{3},t)\,\widetilde{\psi}_{u}(r_{4},t)\;\widetilde{\psi}_{u}^{+}(r_{4},t)\,\widetilde{\psi}_{u}^{+}(r_{3},t)\,\widetilde{\psi}_{d}^{+}(r_{2},t)\,\widetilde{\psi}_{d}^{+}(r_{1},t)}$
$\displaystyle=$
$\displaystyle\frac{1}{(\sqrt{V})^{8}}\mathop{\displaystyle\sum}\limits_{k_{1},k_{2},k_{3},k_{4}}\mathop{\displaystyle\sum}\limits_{k_{5},k_{6},k_{7},k_{8}}\exp
i\\{(k_{1}-k_{8})\cdot r_{1}\\}\;\exp i\\{(k_{2}-k_{7})\cdot r_{2}\\}$
$\displaystyle\times\exp i\\{(k_{3}-k_{6})\cdot r_{3}\\}\;\exp
i\\{(k_{4}-k_{5})\cdot r_{4}\\}$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(k_{1},t)\,\widetilde{\phi}_{d}(k_{2},t)\,\widetilde{\phi}_{u}(k_{3},t)\,\widetilde{\phi}_{u}(k_{4},t)\;\widetilde{\phi}_{u}^{+}(k_{5},t)\,\widetilde{\phi}_{u}^{+}(k_{6},t)\,\widetilde{\phi}_{d}^{+}(k_{7},t)\,\widetilde{\phi}_{d}^{+}(k_{8},t)}$

The case of temperature evolution leads to the same forms with $t$ replaced by
$\beta$ in Eq.(LABEL:Eq.QCFTwoCooperTimeEvn) and in Eq.
(LABEL:Eq.QCFTwoCooperMtmTimeEvn), along with the additional factor $1/Z$ on
the right side of Eq. (LABEL:Eq.QCFTwoCooperMtmTimeEvn) to allow for
$\widehat{\rho}$ being replaced by $\widehat{\sigma}$ when temperature
evolution is involved.

The momentum stochastic correlation function of interest is

$\displaystyle
X(d\,k_{1},d\,k_{2},u\,k_{3},u\,k_{4},u^{+}k_{5},u^{+}k_{6},d^{+}k_{7},d^{+}k_{8})$
(128) $\displaystyle=$
$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{d}(k_{2})\,\widetilde{\phi}_{u}(k_{3})\,\widetilde{\phi}_{u}(k_{4})\;\widetilde{\phi}_{u}^{+}(k_{5})\,\widetilde{\phi}_{u}^{+}(k_{6})\,\widetilde{\phi}_{d}^{+}(k_{7})\,\widetilde{\phi}_{d}^{+}(k_{8})}$

where we have left the $t$ or $\beta$ dependences implicit to shorten the
notation.

###  4.3 Number of Fermions - Partition Function Factor

The mean number of fermions can be obtained by taking the mean of the number
operator (3). The mean equals $N$, the fixed number of fermions. By expanding
the position field operators in terms of plane waves and carrying out the
space integrals we find that for temperature evolution

$\displaystyle\left\langle\widehat{N}\right\rangle$ $\displaystyle=$
$\displaystyle N$ (129) $\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{\alpha}\int
dr\mathbf{\,}\mbox{Tr}\hat{\Psi}_{\alpha}(r)\,\widehat{\sigma}\,\hat{\Psi}_{\alpha}(r)^{{\dagger}}$
$\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{\alpha}\sum_{k}\overline{\widetilde{\phi}_{\alpha}(k,\beta)\,\widetilde{\phi}_{\alpha}^{+}(k,\beta)}$

showing there is a relationship between the partition function $Z$, the number
of fermions $N$ and QCF of the form $X(\alpha k,\alpha
k)=\overline{\widetilde{\phi}_{\alpha}(k,\beta)\,\widetilde{\phi}_{\alpha}^{+}(k\mathbf{,}\beta)}$.
In applying the theory to temperature evolution we have

$\frac{1}{Z}=\frac{N}{\left(\sum_{\alpha}\sum_{k}X(\alpha k,\alpha k)\right)}$
(130)

so the $1/Z$ factors in the results for $X(d\,r_{1},u\,r_{2})$ and
$X(d\,r_{1},d\,r_{2,}\,u\,r_{3},u\,r_{4})$ may be replaced by the last
expression. Hence for example Eq.(124) gives

$\displaystyle\frac{X(d\,r_{1},u\,r_{2};\beta)}{N}$ (131) $\displaystyle=$
$\displaystyle\left\\{\frac{1}{(\sqrt{V})^{4}}\mathop{\displaystyle\sum}\limits_{k_{1},k_{2},k_{3},k_{4}}\exp
i\\{(k_{1}-k_{4})\cdot r_{1}\\}\;\exp i\\{(k_{2}-k_{3})\cdot
r_{2}\\}\times\overline{\widetilde{\phi}_{d}(k_{1},\beta)\,\widetilde{\phi}_{u}(k_{2},\beta)\;\widetilde{\phi}_{u}^{+}(k_{3},\beta)\,\widetilde{\phi}_{d}^{+}(k_{4},\beta)}\right\\}$
$\displaystyle{\Large\div}\left\\{\sum_{\alpha}\sum_{k}\overline{\widetilde{\phi}_{\alpha}(k,\beta)\,\widetilde{\phi}_{\alpha}^{+}(k,\beta)}\right\\}$

so only the stochastic momentum QCF are involved.

##  5 Single Cooper Pair - First Order Changes to QCF

In this Section we will derive expressions for the first order changes in the
QCF for a single Cooper pair based on stochastic momentum fields by carrying
out the stochastic averages analytically The results will be correct to order
$\delta t$ or $\delta\beta$ for time or temperature evolution respectively.

###  5.1 One Cooper Pair - Time Evolution

The first order change in the QCF
$\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
due to time evolution between $t$ and $t+\delta t$ is derived in Appendix 11
and the result is

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1},t+\delta
t)\,\widetilde{\phi}_{u}(k_{2},t+\delta
t)\;\widetilde{\phi}_{u}^{+}(k_{3},t+\delta
t)\,\widetilde{\phi}_{d}^{+}(k_{4},t+\delta t)}$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}-\frac{\hbar^{2}k_{3}^{2}}{2m}-\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$
$\displaystyle+2\lambda^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}\\}\,\delta
t\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$
$\displaystyle-2\lambda^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{(k_{1}+k_{2}),(l_{1}+l_{2})}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\\}\,\delta
t\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$

The key feature of this last result is that the
$\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
at time $t+\delta t$ depend linearly on the
$\overline{\widetilde{\phi}_{d}(l_{1})\,\widetilde{\phi}_{u}(l_{2})\;\widetilde{\phi}_{u}^{+}(l_{3})\,\widetilde{\phi}_{d}^{+}(l_{4})}$
at time $t$, and the change in these QCF is proportional to $\delta t$.
Furthermore, the QCF are c-numbers so their time dependence can be calculated
on a computer. The Grassmann variables are no longer present. The behaviour
depends of course on the initial conditions for
$\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
at $t=0$, and this will be discussed below.

###  5.2 Initial Condition - Time Evolution from t=0

A natural initial condition at $t=0$ would be to assume the fermionic atoms
were non-interacting, and to turn on the interaction via Feshbach resonance
methods. The choice of the interaction constant could span the range for the
BEC/BCS crossover. If the $N$ fermionic atoms are non-interacting the initial
state could be that for zero temperature, in which case all plane wave modes
up to the Fermi surface would be occupied by one spin down atom and one spin
up atom. The density operator is based on a pure state
$\left|\Psi\right\rangle$

$\displaystyle\widehat{\rho}$ $\displaystyle=$
$\displaystyle\left|\Psi\right\rangle\left\langle\Psi\right|$
$\displaystyle\left|\Psi\right\rangle$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\prod}\limits_{k,(|k|\leq
k_{F})}\left|1_{ku}\right\rangle\otimes\mathop{\displaystyle\prod}\limits_{k,(|k|\leq
k_{F})}\left|1_{kd}\right\rangle$ (133)

where $k_{F}=(3\pi^{2}N/V)^{1/3\text{ \thinspace}}$is the radius of the Fermi
sphere.

Hence

$\displaystyle\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,\left|\Psi\right\rangle$
(134) $\displaystyle=$ $\displaystyle\mathop{\displaystyle\prod}\limits_{k\neq
k_{1},(|k|\leq
k_{F})}\left|1_{ku}\right\rangle\,\left|0_{k_{1}u}\right\rangle\otimes\mathop{\displaystyle\prod}\limits_{k\neq
k_{2},(|k|\leq
k_{F})}\left|1_{kd}\right\rangle\,\left|0_{k_{2}u}\right\rangle\qquad
k_{1,}k_{2}\quad inside\;FS$ $\displaystyle=$ $\displaystyle 0\qquad
k_{1,}k_{2}\quad outside\;FS$

so that for $t=0$

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
(135) $\displaystyle=$
$\displaystyle\mbox{Tr}(\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,\hat{\rho}\,\widehat{\Phi}_{u}(k_{3})^{{\dagger}}\,\widehat{\Phi}_{d}(k_{4})^{{\dagger}}).$
$\displaystyle=$
$\displaystyle\delta_{k_{1},k_{4}}\;\delta_{k_{2},k_{3}}\qquad
k_{1,}k_{2},k_{3},k_{4}\quad inside\;FS$ $\displaystyle=$ $\displaystyle
0\qquad otherwise$

Thus from Eq.(LABEL:Eq.QCFOnePairTimeEvnFirstOrderChange) the non-zero
$\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
at $t=\delta t$ that are present in first order must originate from
$\overline{\widetilde{\phi}_{d}(l_{1},0)\,\widetilde{\phi}_{u}(l_{2},0)\;\widetilde{\phi}_{u}^{+}(l_{2},0)\,\widetilde{\phi}_{d}^{+}(l_{1},0)}=1$
with $l_{1},l_{2}$ not outside the Fermi sphere. These are
$\overline{\widetilde{\phi}_{d}(l_{1},\delta
t)\,\widetilde{\phi}_{u}(l_{2},\delta
t)\;\widetilde{\phi}_{u}^{+}(k_{3},\delta
t)\,\widetilde{\phi}_{d}^{+}(k_{4},\delta t)}$ where $k_{3}+k_{4}=l_{2}+l_{1}$
from the first $2\lambda^{2}$ term and
$\overline{\widetilde{\phi}_{d}(k_{1},\delta
t)\,\widetilde{\phi}_{u}(k_{2},\delta
t)\;\widetilde{\phi}_{u}^{+}(l_{2},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{1},\delta t)}$ where $k_{1}+k_{2}=l_{1}+l_{2}$
from the second $2\lambda^{2}$ term These restrictions represent momentum
conservation.

A specific case of interest is where $l_{1}=l_{4}=+k$ and $l_{2}=l_{3}=-k$,
where $k$ is on the Fermi sphere. The term
$\overline{\widetilde{\phi}_{d}(+k,\delta t)\,\widetilde{\phi}_{u}(-k,\delta
t)\;\widetilde{\phi}_{u}^{+}(+k,\delta t)\,\widetilde{\phi}_{d}^{+}(-k,\delta
t)}$ could originate from the first $2\lambda^{2}\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term but would not originate from the second $2\lambda^{2}$ $\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term Hence there would be no cancellation of the contributions. On the other
hand, for the case where $l_{1}=+k$, $l_{2}=-k$, $l_{3}=-k$ and $l_{4}=+k$
where $k$ is on the Fermi sphere, the term
$\overline{\widetilde{\phi}_{d}(+k,\delta t)\,\widetilde{\phi}_{u}(-k,\delta
t)\;\widetilde{\phi}_{u}^{+}(-k,\delta t)\,\widetilde{\phi}_{d}^{+}(+k,\delta
t)}$ could originate from the first $2\lambda^{2}\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term and could also originate from the second $2\lambda^{2}$ $\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term. However in this case the two contributions cancell out.

###  5.3 One Cooper Pair - Temperature Evolution

The first order change in the QCF
$\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
due to temperature evolution between $\beta$ and $\beta+\delta\beta$ is
derived in Appendix 11 and the result is

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1},\beta+\delta\beta)\,\widetilde{\phi}_{u}(k_{2},\beta+\delta\beta)\;\widetilde{\phi}_{u}^{+}(k_{3},\beta+\delta\beta)\,\widetilde{\phi}_{d}^{+}(k_{4},\beta+\delta\beta)}$
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}+\frac{\hbar^{2}k_{3}^{2}}{2m}+\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$
$\displaystyle-2\eta^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}\\}\,\delta\beta\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$
$\displaystyle-2\eta^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{(k_{1}+k_{2}),(l_{1}+l_{2})}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\\}\,\delta\beta\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$

Similar comments to the time evolution case apply here as well. The initial
conditions are discussed below.

###  5.4 Initial Condition - Temperature Evolution from $\beta=0$

A natural initial condition at $\beta=0$ would be to assume the fermionic
atoms were non-interacting, and to turn on the interaction via Feshbach
resonance methods. The choice of the interaction constant could span the range
for the BEC/BCS crossover. If the $N$ fermionic atoms are non-interacting the
initial state could be that for a very high temperature, in which case a large
number of plane wave modes could be occupied, both inside and outside the
Fermi surface by at most one spin down atom and one spin up atom. The density
operator is based on a mixed state. Suppose we consider the case where there
are $n$ different modes $k$, where $|k|$ would range from being $\ll k_{F}$ to
being $\gg k_{F}$. The one-fermion states are $\left|1_{ku}\right\rangle$ or
$\left|1_{kd}\right\rangle$ if occupied, and $\left|0_{ku}\right\rangle$ or
$\left|0_{kd}\right\rangle$ if unoccupied. So a typical $N$ fermion state will
be
$\left|\lambda\right\rangle=\mathop{\displaystyle\prod}\limits_{k}\;\left|\mu_{ku}\right\rangle\left|\nu_{kd}\right\rangle$
where $\mu,\nu=0,1$ only and with $\lambda\equiv\\{\mu_{ku},\nu_{kd}\\}$. The
number of distinct states $\left|\lambda\right\rangle$ is given by the number
of ways of choosing $N$ occupied one fermion states out of the set of $2n$
possible one fermion states, and this is given by
$\mathcal{N=}^{2n}C_{N}=2n!/[(N)!(2n-N)!]\approx(2n)^{N}/N!.$since we have
$n\gg N$. At high temperatures each of these states
$\left|\lambda\right\rangle$ is equally probable, so the density operator is
given by

$\displaystyle\widehat{\rho}$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\mathop{\displaystyle\sum}\limits_{\lambda}\left|\lambda\right\rangle\left\langle\lambda\right|$
(137) $\displaystyle\approx$
$\displaystyle\frac{N!}{(2n)^{N}}\mathop{\displaystyle\sum}\limits_{\lambda}\left|\lambda\right\rangle\left\langle\lambda\right|$

We then want to consider the quantity

$\displaystyle\mbox{Tr}(\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,\hat{\rho}\,\widehat{\Phi}_{u}(k_{3})^{{\dagger}}\,\widehat{\Phi}_{d}(k_{4})^{{\dagger}})$
(138) $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\mathop{\displaystyle\sum}\limits_{\lambda}\mbox{Tr}(\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,(\left|\lambda\right\rangle\left\langle\lambda\right|)\,\widehat{\Phi}_{u}(k_{3})^{{\dagger}}\,\widehat{\Phi}_{d}(k_{4})^{{\dagger}})$

so to get a non-zero contribution for a given $\lambda$ we need
$\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,\left|\lambda\right\rangle\neq
0$ and
$\widehat{\Phi}_{d}(k_{4})\,\widehat{\Phi}_{u}(k_{3})\,\left|\lambda\right\rangle\neq
0$. Noting that $\widehat{\Phi}_{\alpha}(k)\left|1_{k\alpha}\right\rangle$
$=\left|0_{k\alpha}\right\rangle\,$and
$\widehat{\Phi}_{\alpha}(k)\left|0_{k\alpha}\right\rangle$ $=0$ we see that to
get a no-zero contribution for a given $\lambda$ we require $k_{1}=k_{4}$ and
$k_{2}=k_{3}$ which restricts $\lambda$ to be of the form
$\left|\lambda^{\\#}\right\rangle=\left|1_{k_{1}d}\right\rangle\left|1_{k_{2}u}\right\rangle\mathop{\displaystyle\prod}\limits_{k\neq
k_{1,}k_{2}}\;\left|\mu_{ku}\right\rangle\left|\nu_{kd}\right\rangle$.
Similarly to the previous situation, there are
$\mathcal{N}^{\\#}=^{2n-2}C_{N-2}\approx(2n)^{(N-2)}/(N-2)!$ different states
$\left|\lambda^{\\#}\right\rangle$ (since this equals the number of ways of
choosing $N-2$ occupied one fermion states from $2n-2$ possible remaining one
fetmi0n states). In this case we have

$\displaystyle\mbox{Tr}(\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,\hat{\rho}\,\widehat{\Phi}_{u}(k_{3})^{{\dagger}}\,\widehat{\Phi}_{d}(k_{4})^{{\dagger}})$
(139) $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\delta_{k_{1,}k_{4}}\,\delta_{k_{2,}k_{3}}\,\mathop{\displaystyle\sum}\limits_{\lambda^{\\#}}\mbox{Tr}(\,\left|0_{k_{1}d}\right\rangle\left|0_{k_{2}u}\right\rangle\left(\mathop{\displaystyle\prod}\limits_{k\neq
k_{1,}k_{2}}\;\left|\mu_{ku}\right\rangle\left|\nu_{kd}\right\rangle\otimes\left\langle\mu_{ku}\right|\left\langle\nu_{kd}\right|\right)\,\left\langle
0_{k_{1}d}\right|\left\langle 0_{k_{2}u}\right|\,)$ $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{N}}\delta_{k_{1,}k_{4}}\,\delta_{k_{2,}k_{3}}\,\mathop{\displaystyle\sum}\limits_{\lambda^{\\#}}\,1=\frac{\mathcal{N}^{\\#}}{\mathcal{N}}\delta_{k_{1,}k_{4}}\,\delta_{k_{2,}k_{3}}\,$
$\displaystyle\approx$
$\displaystyle\frac{1}{4}\left(\frac{N}{n}\right)^{2}\delta_{k_{1,}k_{4}}\,\delta_{k_{2,}k_{3}}$

Writing $\widehat{\sigma}=Z\,\widehat{\rho}$ we then see that at $\beta=0$

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
(140) $\displaystyle=$
$\displaystyle\mbox{Tr}(\widehat{\Phi}_{d}(k_{1})\,\widehat{\Phi}_{u}(k_{2})\,\widehat{\sigma}\,\widehat{\Phi}_{u}(k_{3})^{{\dagger}}\,\widehat{\Phi}_{d}(k_{4})^{{\dagger}}).$
$\displaystyle\approx$
$\displaystyle\frac{Z}{4}\left(\frac{N}{n}\right)^{2}\delta_{k_{1,}k_{4}}\,\delta_{k_{2,}k_{3}}$

Thus from Eq.(LABEL:Eq.QCFOnePairTempEvnFirstOrderChangeA) the non-zero
$\overline{\widetilde{\phi}_{d}(k_{1})\,\widetilde{\phi}_{u}(k_{2})\;\widetilde{\phi}_{u}^{+}(k_{3})\,\widetilde{\phi}_{d}^{+}(k_{4})}$
at $\beta=\delta\beta$ that are present in first order must originate from
$\overline{\widetilde{\phi}_{d}(l_{1},0)\,\widetilde{\phi}_{u}(l_{2},0)\;\widetilde{\phi}_{u}^{+}(l_{2},0)\,\widetilde{\phi}_{d}^{+}(l_{1},0)}\approx\frac{Z}{4}\left(\frac{N}{n}\right)^{2}$
with $l_{1},l_{2}$ not restricted to the Fermi sphere. These are
$\overline{\widetilde{\phi}_{d}(l_{1},\delta\beta)\,\widetilde{\phi}_{u}(l_{2},\delta\beta)\;\widetilde{\phi}_{u}^{+}(k_{3},\delta\beta)\,\widetilde{\phi}_{d}^{+}(k_{4},\delta\beta)}$
where $k_{3}+k_{4}=l_{2}+l_{1}$ from the first $2\eta^{2}$ term and
$\overline{\widetilde{\phi}_{d}(k_{1},\delta\beta)\,\widetilde{\phi}_{u}(k_{2},\delta\beta)\;\widetilde{\phi}_{u}^{+}(l_{2},\delta\beta)\,\widetilde{\phi}_{d}^{+}(l_{1},\delta\beta)}$
where $k_{1}+k_{2}=l_{1}+l_{2}$ from the second $2\eta^{2}$ term These
restrictions represent momentum conservation.

A specific case of interest is where $l_{1}=l_{4}=+k$ and $l_{2}=l_{3}=-k$,
where $k$ need not be on the Fermi sphere. The term
$\overline{\widetilde{\phi}_{d}(+k,\delta\beta)\,\widetilde{\phi}_{u}(-k,\delta\beta)\;\widetilde{\phi}_{u}^{+}(+k,\delta\beta)\,\widetilde{\phi}_{d}^{+}(-k,\delta\beta)}$
could originate from the first $2\eta^{2}\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term but would not originate from the second $2\eta^{2}$ $\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term Hence there would be no cancellation of the contributions. On the other
hand, for the case where $l_{1}=+k$, $l_{2}=-k$, $l_{3}=-k$ and $l_{4}=+k$
where $k$ need not be on the Fermi sphere, the term
$\overline{\widetilde{\phi}_{d}(+k,\delta\beta)\,\widetilde{\phi}_{u}(-k,\delta\beta)\;\widetilde{\phi}_{u}^{+}(-k,\delta\beta)\,\widetilde{\phi}_{d}^{+}(+k,\delta\beta)}$
could originate from the first $2\eta^{2}\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term and could also originate from the second $2\eta^{2}$ $\times$
$\overline{\widetilde{\phi}_{d}(+k,0)\,\widetilde{\phi}_{u}(-k,0)\;\widetilde{\phi}_{u}^{+}(-k,0)\,\widetilde{\phi}_{d}^{+}(+k,0)}$
term. However in this temperature evolution case the two contributions do not
cancell out - in contrast to the time evolution case.

##  6 Two Cooper Pairs - First Order Changes to QCF

In this Section we will derive expressions for the first order changes in the
QCF for two Cooper pairs based on stochastic momentum fields by carrying out
the stochastic averages analytically The results will be correct to order
$\delta t$ or $\delta\beta$ for time or temperature evolution respectively.

###  6.1 Two Cooper Pairs - Time Evolution

For the case of time evolution the QCF for two Cooper pairs at time $t+\delta
t$ is related to that at time $t$ via

$\displaystyle\left[\begin{array}[]{c}\widetilde{\phi}_{d}(k_{1},t+\delta
t)\,\widetilde{\phi}_{d}(k_{2},t+\delta
t)\;\widetilde{\phi}_{u}(k_{3},t+\delta
t)\,\widetilde{\phi}_{u}(k_{4},t+\delta t)\\\
\times\;\widetilde{\phi}_{u}^{+}(k_{5},t+\delta
t)\,\widetilde{\phi}_{u}^{+}(k_{6},t+\delta
t)\widetilde{\phi}_{d}^{+}(k_{7},t+\delta
t)\,\widetilde{\phi}_{d}^{+}(k_{8},t+\delta t)\end{array}\right]_{StochAver}$
(159) $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}}\mathop{\displaystyle\sum}\limits_{l_{2}}\mathop{\displaystyle\sum}\limits_{l_{3}}\mathop{\displaystyle\sum}\limits_{l_{4}}\mathop{\displaystyle\sum}\limits_{l_{5}}\mathop{\displaystyle\sum}\limits_{l_{6}}\mathop{\displaystyle\sum}\limits_{l_{7}}\mathop{\displaystyle\sum}\limits_{l_{8}}$
$\displaystyle\left[\begin{array}[]{c}(\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\,\delta_{k_{5},l_{5}}\,\delta_{k_{6},l_{6}}\,\delta_{k_{7},l_{7}}\,\delta_{k_{8},l_{8}})\\\
\times\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}+\frac{\hbar^{2}k_{3}^{2}}{2m}+\frac{\hbar^{2}k_{4}^{2}}{2m}-\frac{\hbar^{2}k_{5}^{2}}{2m}-\frac{\hbar^{2}k_{6}^{2}}{2m}-\frac{\hbar^{2}k_{7}^{2}}{2m}-\frac{\hbar^{2}k_{8}^{2}}{2m}\right\\}\delta
t\right)\\\ \\\
+(-2\lambda^{2})\,\delta_{k_{1},l_{1}}\delta_{k_{2},l_{2}}\delta_{k_{3},l_{3}}\delta_{k_{4},l_{4}}\,\delta
t\\\
\times\left(\begin{array}[]{c}\delta_{k_{6},l_{5}}\delta_{k_{8},l_{8}}\;\delta_{(k_{5}+k_{7}),(l_{6}+l_{7})}+\delta_{k_{6},l_{5}}\delta_{k_{7},l_{7}}\;\delta_{(k_{5}+k_{8}),(l_{6}+l_{8})}\\\
+\delta_{k_{5},l_{6}}\delta_{k_{8},l_{8}}\;\delta_{(k_{6}+k_{7}),(l_{5}+l_{7})}+\delta_{k_{5},l_{6}}\delta_{k_{7},l_{7}}\;\delta_{(k_{6}+k_{8}),(l_{5}+l_{8})}\end{array}\right)\\\
\\\
+(2\lambda^{2})\,\delta_{k_{5},l_{5}}\,\delta_{k_{6},l_{6}}\,\delta_{k_{7},l_{7}}\,\delta_{k_{8},l_{8}}\,\delta
t\\\
\times\left(\begin{array}[]{c}\delta_{k_{1},l_{1}}\delta_{k_{3},l_{4}}\;\delta_{(k_{2}+k_{4}),(l_{2}+l_{3})}+\delta_{k_{1},l_{1}}\delta_{k_{4},l_{3}}\;\delta_{(k_{2}+k_{3}),(l_{2}+l_{4})}\\\
+\delta_{k_{2},l_{2}}\delta_{k_{3},l_{4}}\;\delta_{(k_{1}+k_{4}),(l_{1}+l_{3})}+\delta_{k_{2},l_{2}}\delta_{k_{4},l_{3}}\;\delta_{(k_{1}+k_{3}),(l_{1}+l_{4})}\end{array}\right)\end{array}\right]$
$\displaystyle\times\left[\begin{array}[]{c}\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{d}(l_{2},t)\;\widetilde{\phi}_{u}(l_{3},t)\,\widetilde{\phi}_{u}(l_{4},t)\\\
\times\;\widetilde{\phi}_{u}^{+}(l_{5},t)\,\widetilde{\phi}_{u}^{+}(l_{6},t)\widetilde{\phi}_{d}^{+}(l_{7},t)\,\widetilde{\phi}_{d}^{+}(l_{8},t)\end{array}\right]_{StochAver}$

Thus the first order change in the QCF for two Cooper pairs depends linearly
on $\delta t$. The derivation is given in Sect.11.3 Note the signs for the
$2\eta^{2}$ terms are opposite and that halk the kinetetic energy terms
$\frac{\hbar^{2}k^{2}}{2m}$ have opposite signs. Similar intial conditions to
the single Cooper pair time evolution apply.

###  6.2 Two Cooper Pairs - Temperature Evolution

For the case of temperature evolution the QCF for two Cooper pairs at
temperature $\beta+\delta\beta$ is related to that at temperature $\beta$ via

$\displaystyle\left[\begin{array}[]{c}\widetilde{\phi}_{d}(k_{1},\beta+\delta\beta)\,\widetilde{\phi}_{d}(k_{2},\beta+\delta\beta)\;\widetilde{\phi}_{u}(k_{3},\beta+\delta\beta)\,\widetilde{\phi}_{u}(k_{4},\beta+\delta\beta)\\\
\times\;\widetilde{\phi}_{u}^{+}(k_{5},\beta+\delta\beta)\,\widetilde{\phi}_{u}^{+}(k_{6},\beta+\delta\beta)\widetilde{\phi}_{d}^{+}(k_{7},\beta+\delta\beta)\,\widetilde{\phi}_{d}^{+}(k_{8},\beta+\delta\beta)\end{array}\right]_{StochAver}$
(179) $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}}\mathop{\displaystyle\sum}\limits_{l_{2}}\mathop{\displaystyle\sum}\limits_{l_{3}}\mathop{\displaystyle\sum}\limits_{l_{4}}\mathop{\displaystyle\sum}\limits_{l_{5}}\mathop{\displaystyle\sum}\limits_{l_{6}}\mathop{\displaystyle\sum}\limits_{l_{7}}\mathop{\displaystyle\sum}\limits_{l_{8}}$
$\displaystyle\left[\begin{array}[]{c}(\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\,\delta_{k_{5},l_{5}}\,\delta_{k_{6},l_{6}}\,\delta_{k_{7},l_{7}}\,\delta_{k_{8},l_{8}})\\\
\times\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}+\frac{\hbar^{2}k_{3}^{2}}{2m}+\frac{\hbar^{2}k_{4}^{2}}{2m}+\frac{\hbar^{2}k_{5}^{2}}{2m}+\frac{\hbar^{2}k_{6}^{2}}{2m}+\frac{\hbar^{2}k_{7}^{2}}{2m}+\frac{\hbar^{2}k_{8}^{2}}{2m}\right\\}\delta\beta\right)\\\
\\\
+(+2\eta^{2})\,\delta_{k_{1},l_{1}}\delta_{k_{2},l_{2}}\delta_{k_{3},l_{3}}\delta_{k_{4},l_{4}}\,\delta\beta\\\
\times\left(\begin{array}[]{c}\delta_{k_{6},l_{5}}\delta_{k_{8},l_{8}}\;\delta_{(k_{5}+k_{7}),(l_{6}+l_{7})}+\delta_{k_{6},l_{5}}\delta_{k_{7},l_{7}}\;\delta_{(k_{5}+k_{8}),(l_{6}+l_{8})}\\\
+\delta_{k_{5},l_{6}}\delta_{k_{8},l_{8}}\;\delta_{(k_{6}+k_{7}),(l_{5}+l_{7})}+\delta_{k_{5},l_{6}}\delta_{k_{7},l_{7}}\;\delta_{(k_{6}+k_{8}),(l_{5}+l_{8})}\end{array}\right)\\\
\\\
+(2\eta^{2})\,\delta_{k_{5},l_{5}}\,\delta_{k_{6},l_{6}}\,\delta_{k_{7},l_{7}}\,\delta_{k_{8},l_{8}}\,\delta\beta\\\
\times\left(\begin{array}[]{c}\delta_{k_{1},l_{1}}\delta_{k_{3},l_{4}}\;\delta_{(k_{2}+k_{4}),(l_{2}+l_{3})}+\delta_{k_{1},l_{1}}\delta_{k_{4},l_{3}}\;\delta_{(k_{2}+k_{3}),(l_{2}+l_{4})}\\\
+\delta_{k_{2},l_{2}}\delta_{k_{3},l_{4}}\;\delta_{(k_{1}+k_{4}),(l_{1}+l_{3})}+\delta_{k_{2},l_{2}}\delta_{k_{4},l_{3}}\;\delta_{(k_{1}+k_{3}),(l_{1}+l_{4})}\end{array}\right)\end{array}\right]$
$\displaystyle\times\left[\begin{array}[]{c}\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{d}(l_{2},\beta)\;\widetilde{\phi}_{u}(l_{3},\beta)\,\widetilde{\phi}_{u}(l_{4},\beta)\\\
\times\;\widetilde{\phi}_{u}^{+}(l_{5},\beta)\,\widetilde{\phi}_{u}^{+}(l_{6},\beta)\widetilde{\phi}_{d}^{+}(l_{7},\beta)\,\widetilde{\phi}_{d}^{+}(l_{8},\beta)\end{array}\right]_{StochAver}$

Thus the first order change in the QCF for two Cooper pairs depends linearly
on $\delta\beta$. The derivation is similar to the case of time evolution.
Note the signs for the $2\eta^{2}$ terms are now the same and that all the
kinetetic energy terms $\frac{\hbar^{2}k^{2}}{2m}$ also have the same sign.
Similar intial conditions to the single Cooper pair temperature evolution also
apply.

##  7 Solution of Evolution Equations for QCF

The evolution equations for the QCF can be written in terms of the definitions
in Eqs.(125) and (128) for the one Cooper pair QCF
$X(d\,k_{1},u\,k_{2},u^{+}\,k_{3},d^{+}\,k_{4})$ and the two Cooper pair QCF
$X(d\,k_{1},d\,k_{2},u\,k_{3},u\,k_{4},u^{+}k_{5},u^{+}k_{6},d^{+}k_{7},d^{+}k_{8})$
in the forms

$\displaystyle\frac{\partial}{\partial\xi}X(d\,k_{1},u\,k_{2},u^{+}\,k_{3},d^{+}\,k_{4})$
(182) $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}M(k_{1},k_{2},k_{3},k_{4};l_{1},l_{2},l_{3},l_{4})\,X(d\,l_{1},u\,l_{2},u^{+}\,l_{3},d^{+}\,l_{4})$
$\displaystyle\frac{\partial}{\partial\xi}X(d\,k_{1},d\,k_{2},u\,k_{3},u\,k_{4},u^{+}k_{5},u^{+}k_{6},d^{+}k_{7},d^{+}k_{8})$
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}M(k_{1},k_{2},k_{3},k_{4},k_{5},k_{6},k_{7},k_{8};l_{1},l_{2},l_{3},l_{4},l_{5},l_{6},l_{7,}l_{8})\,$
$\displaystyle\times
X(d\,l_{1},d\,l_{2},u\,l_{3},u\,l_{4},u^{+}l_{5},u^{+}l_{6},d^{+}l_{7},d^{+}l_{8})$

where $\xi$ is either $t$ or $\beta$ for time or temperature evolution, and
the matrices $M$ can be read off from Eq.
(LABEL:Eq.QCFOnePairTimeEvnFirstOrderChange),
(LABEL:Eq.QCFOnePairTempEvnFirstOrderChangeA) or
(LABEL:Eq.QCFOTwoPairsTempEvnFirstOrderChangeA). By inspection it can be seen
that $M$ is a real, symmetrric matrix.

These equations can be solved using the normalised, orthogonal column
eigenvectors $x_{i}$ of $M$. Thus with real eigenvalues $\lambda_{i}$

$M\,x_{i}=\lambda_{i}\,x_{i}\qquad x_{i}^{T}\;x_{j}=\delta_{ij}$ (183)

the solution is of the form

$\displaystyle X(\xi)$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{i}c_{i}(0)\;\exp(\lambda_{i}\,\xi)\;x_{i}$
(184) $\displaystyle c_{i}(0)$ $\displaystyle=$ $\displaystyle
x_{i}^{T}\;X(0)$ (185)

where the initial condition is specified by the coefficients $c_{i}(0)$ which
are determined from the initial column vector $X(0)$ .

This method of solution requires first calculating the matrix $M$ and the
initial column vector $X(0)$, then determining the eigenvectors and
eigenvalues of $M$. The alternative is a step-by-step development of the QCF
using Eq.(LABEL:Eq.QCFOnePairTimeEvnFirstOrderChange),
(LABEL:Eq.QCFOnePairTempEvnFirstOrderChangeA) or
(LABEL:Eq.QCFOTwoPairsTempEvnFirstOrderChangeA)

##  8 Summary

In the present paper it has been shown how Grassmann Phase Space Theory (GSPT)
can be applied to the BEC/BCS crossover in cold fermionic atomic gases and
used to determine the evolution (over either time or temperature) of the
Quantum Correlation Functions (QCF) that specify: (a) the positions of the
spin up and spin down fermionic atoms in a single Cooper pair and (b) the
positions of the two spin up and two spin down fermionic atoms in two Cooper
pairs The first of these QCF is relevent to describing the change in size of a
Cooper pair, as the fermion-fermion coupling constant is changed through the
crossover from a small Cooper pair on the BEC side to a large Cooper pair on
the BCS side. The second of these QCF is important for describing the
correlations berween the positions of the fermionic atoms in two Cooper pairs,
which is expected to be small at the BEC or BCS sides of the crossover, but is
expected to be significant in the strong interaction unitary regime, where the
size of a Cooper pair is comparable to the separation between Cooper pairs. In
GPST the QCF are given as the stochastic average of products of Grassman
stochastic position fields, which are related via Fourier transforms to the
stochastic average of products of Grassman stochastic momentum fields. Using
the no-correlation theorem, GPST shows that the stochastic average of the
products of Grassman stochastic momentum fields at a later time (or lower
temperature) is related linearly to the stochastic average of the products of
Grassman stochastic momentum fields at an earlier time (or higher
temperature), and furthermore that the matrix elements involved in the linear
relations are all c-numbers. Expressions for these matrix elements
corresponding to a small time increment (or a small temperature change) have
been obtained analytically, providing the formulae needed for numerical
studies of the evolution that are planned for a future publication. The
initial conditions envisaged include those for a non-interacting fermionic gas
at zero temperature to study the time evolution of the creation of a Cooper
pair when the interaction is switched on via Feshbach resonance methods. Other
initial conditions described include a high temperature gas (where the effect
of the interactions can be ignored in the initial state), and where the
evolution will be studied as the temperature is lowered until either a BEC
state with small Cooper pairs or a BCS state with large Cooper pairs forms,
depending on the fermion-fermion coupling constant. The behaviour for the case
where the coupling constant is very large (the unitary regime) would be of
particular interest. Full derivations of the expressions have been presented
in the Appendices.

##  9 Acknowledgements

This research was supported by the Centre for Quantum Technology Theory,
School of Science, Computing and Engineering Technology, Swinburne University,
Melbourne, Victoria 3122, Australia.

##  10 Appendix A - Details for Grassmann Field Theory

In this Appendix some details for the Grassmann field theory Section are
provided.

###  10.1 Gaussian-Markoff Noise Terms

The basic stochastic average properties of the Gaussian-Markoff random noise
terms are

$\displaystyle\overline{\Gamma_{{\small a}}(t_{1})}\,$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle\overline{\Gamma_{{\small
a}}(t_{1})\Gamma_{b}(t_{2})}\,$ $\displaystyle=$ $\displaystyle\delta_{{\small
ab}}\delta(t_{1}-t_{2})$
$\displaystyle\overline{\Gamma_{a}(t_{1})\Gamma_{b}(t_{2})\Gamma_{c}(t_{3})}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\overline{\Gamma_{a}(t_{1})\Gamma_{b}(t_{2})\Gamma_{c}(t_{3})\Gamma_{d}(t_{4})}$
$\displaystyle=$
$\displaystyle\overline{\Gamma_{a}(t_{1})\Gamma_{b}(t_{2}})\overline{\Gamma_{c}(t_{3})\Gamma_{d}(t_{4}})+\overline{\Gamma_{a}(t_{1})\Gamma_{c}(t_{3}})\overline{\Gamma_{b}(t_{2})\Gamma_{d}(t_{4}})$
$\displaystyle+\overline{\Gamma_{a}(t_{1})\Gamma_{d}(t_{4}})\overline{\Gamma_{b}(t_{2})\Gamma_{c}(t_{3}})$

showing that the stochastic averages of a single $\Gamma$ is zero and the
stochastic average of the product of two $\Gamma$’s is zero if they are
different and delta function correlated in the time difference if they are the
same. In addition, the stochastic averages of products of odd numbers of
$\Gamma$ are zero and stochastic averages of products of even numbers of
$\Gamma$ are the sums of products of stochastic averages of pairs of $\Gamma$.

###  10.2 Determining the B Matrix - Time Evolution Case

Using the notation in Section 3.6 and after substituting for $\delta(r-s)$,
the non-zero elements of the diffusion matrix in Eq. (21) are given by

$\displaystyle D_{u1,d1}(\psi(s),\psi(r))$ $\displaystyle=$
$\displaystyle\frac{ig}{\hbar
V}\sum\limits_{q}\theta_{d1}^{q}(s)\;\theta_{u1}^{q}(r)\qquad
D_{d1,u1}(\psi(s),\psi(r))=\frac{ig}{\hbar
V}\sum\limits_{q}\theta_{u1}^{q}(s)\;\theta_{d1}^{q}(r)$ $\displaystyle
D_{u2,d2}(\psi(s),\psi(r))$ $\displaystyle=$ $\displaystyle\frac{ig}{\hbar
V}\sum\limits_{q}\theta_{d2}^{q}(s)\;\theta_{u2}^{q}(r)\qquad
D_{d2,u2}(\psi(s),\psi(r))=\frac{ig}{\hbar
V}\sum\limits_{q}\theta_{u2}^{q}(s)\;\theta_{d2}^{q}(r)$

where

$\displaystyle\theta_{d1}^{q}(s)$ $\displaystyle=$
$\displaystyle\psi_{d1}(s)\,\exp(+iqs)\qquad\theta_{u1}^{q}(s)=\psi_{u1}(s)\,\exp(-iqs)$
$\displaystyle\theta_{d2}^{q}(s)$ $\displaystyle=$
$\displaystyle\psi_{d2}(s)\,\exp(-iqs)\qquad\theta_{u2}^{q}(s)=\psi_{u1}(s)\,\exp(+iqs)$
(188)

We write the diffusion matrix in the form

$D_{\alpha A\,\beta
B}[\psi(s,t),s;\psi(r,t),r]=\mathop{\displaystyle\sum}\limits_{q}\mathop{\displaystyle\sum}\limits_{\gamma
C,\delta D}M_{\alpha A,\gamma C;\beta B,\delta D}\;\theta_{\gamma
C}^{q}(s)\,\theta_{\delta D}^{q}(r)$ (189)

where the symmetric matrix $M$ is

$[M]=\frac{ig}{\hbar V}\,\left[\begin{tabular}[]{lllll}$(\alpha A,\gamma
C)\downarrow\backslash\;(\beta B,\delta
D)\rightarrow$&$u1,d1$&$d1,u1$&$u2,d2$&$d2,u2$\\\ $u1,d1$&0&1&0&0\\\
$d1,u1$&1&0&0&0\\\ $u2,d2$&0&0&0&-1\\\ $d2,u2$&0&0&-1&0\end{tabular}\right]$
(190)

Using Takagi [31] factorisation $M=KK^{T}$ the matrix $M$ is given by

$M_{\alpha A,\gamma C;\beta B,\delta
D}=\mathop{\displaystyle\sum}\limits_{\xi}K_{\alpha A,\gamma C;\xi}\,K_{\beta
B,\delta D;\xi}$ (191)

where

$[K]=\sqrt{\frac{ig}{\hbar V}}\,\left[\begin{tabular}[]{lllll}$(\alpha
A,\gamma
C)\downarrow\backslash\;(\xi)\rightarrow$&$u1,d1$&$d1,u1$&$u2,d2$&$d2,u2$\\\
$u1,d1$&+$\frac{\text{1}}{\sqrt{\text{2}}}$&+$\frac{\text{i}}{\sqrt{\text{2}}}$&0&0\\\
$d1,u1$&+$\frac{\text{1}}{\sqrt{\text{2}}}$&-$\,\frac{\text{i}}{\sqrt{\text{2}}}$&0&0\\\
$u2,d2$&0&0&+$\frac{\text{1}}{\sqrt{\text{2}}}$&+$\frac{\text{i}}{\sqrt{\text{2}}}$\\\
$d2,u2$&0&0&-$\frac{\text{1}}{\sqrt{\text{2}}}$&+$\frac{\text{i}}{\sqrt{\text{2}}}$\end{tabular}\right]$
(192)

Thus we have

$\displaystyle D_{\alpha A\,\beta B}[\psi(s,t),s;\psi(r,t),r]$ (193)
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{q}\mathop{\displaystyle\sum}\limits_{\xi}\,\left(\mathop{\displaystyle\sum}\limits_{\gamma
C}K_{\alpha A,\gamma C;\xi}\,\theta_{\gamma
C}^{q}(s)\right)\,\left(\mathop{\displaystyle\sum}\limits_{\delta D}K_{\beta
B,\delta D;\xi}\theta_{\delta D}^{q}(r)\right)$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{q,\xi}B_{q,\xi}^{\alpha
A}[\psi(s,t),s]\,B_{q,\xi}^{\beta B}[\psi(r,t),r]$

where

$B_{q,\xi}^{\alpha A}[\psi(s,t),s]\,=\mathop{\displaystyle\sum}\limits_{\gamma
C}K_{\alpha A,\gamma C;\xi}\,\theta_{\gamma C}^{q}(s)$ (194)

Hence we have written the diffusion matrix in the form $D=BB^{T}$ as required
to determine the Ito SFE.

On substituting for $K_{\alpha A,\gamma C;\xi}$ and $\theta_{\gamma C}^{q}(s)$
from Eqs. (192) and (188) we obtain the Ito SFE in the form given in Eq.
(LABEL:Eq.ThetaMatrixTempEvn).

###  10.3 Determining the B Matrix - Temperature Evolution Case

The temperature evolution case is treated similarly to the time evolution case
in the previous Section.

However in this case the diffusion matrix is

$\displaystyle D_{u1,d1}(\psi(s),\psi(r))$ $\displaystyle=$
$\displaystyle\frac{g}{4V}\sum\limits_{q}\theta_{d1}^{q}(s)\;\theta_{u1}^{q}(r)\qquad
D_{d1,u1}(\psi(s),\psi(r))=\frac{g}{4V}\sum\limits_{q}\theta_{u1}^{q}(s)\;\theta_{d1}^{q}(r)$
$\displaystyle D_{u2,d2}(\psi(s),\psi(r))$ $\displaystyle=$
$\displaystyle\frac{g}{4V}\sum\limits_{q}\theta_{d2}^{q}(s)\;\theta_{u2}^{q}(r)\qquad
D_{d2,u2}(\psi(s),\psi(r))=\frac{g}{4V}\sum\limits_{q}\theta_{u2}^{q}(s)\;\theta_{d2}^{q}(r)$

and the $M$ matrix is

$[M]=\frac{g}{4V}\,\left[\begin{tabular}[]{lllll}$(\alpha A,\gamma
C)\downarrow\backslash\;(\beta B,\delta
D)\rightarrow$&$u1,d1$&$d1,u1$&$u2,d2$&$d2,u2$\\\ $u1,d1$&0&1&0&0\\\
$d1,u1$&1&0&0&0\\\ $u2,d2$&0&0&0&1\\\ $d2,u2$&0&0&1&0\end{tabular}\right]$
(196)

which results in the $K$ matrix

$[K]=\sqrt{\frac{g}{4V}}\,\left[\begin{tabular}[]{lllll}$(\alpha A,\gamma
C)\downarrow\backslash\;(\xi)\rightarrow$&$u1,d1$&$d1,u1$&$u2,d2$&$d2,u2$\\\
$u1,d1$&+$\frac{\text{1}}{\sqrt{\text{2}}}$&+$\frac{\text{i}}{\sqrt{\text{2}}}$&0&0\\\
$d1,u1$&+$\frac{\text{1}}{\sqrt{\text{2}}}$&-$\,\frac{\text{i}}{\sqrt{\text{2}}}$&0&0\\\
$u2,d2$&0&0&+$\frac{\text{1}}{\sqrt{\text{2}}}$&+$\frac{\text{i}}{\sqrt{\text{2}}}$\\\
$d2,u2$&0&0&+$\frac{\text{1}}{\sqrt{\text{2}}}$&-$\frac{\text{i}}{\sqrt{\text{2}}}$\end{tabular}\right]$
(197)

On substituting for $K_{\alpha A,\gamma C;\xi}$ and $\theta_{\gamma C}^{q}(s)$
from Eqs. (197) and (188) we obtain the Ito SFE in the form given in Eq.
(LABEL:Eq.ItoSFETimeEvoln).

###  10.4 Stochastic Momentum Fields

The inverse transformation between the stochastic position and momentum fields
is

$\displaystyle\tilde{\phi}_{\alpha}(k)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{V}}\int ds\;\exp(-ik\cdot
s)\;\tilde{\psi}_{\alpha}(s)$ $\displaystyle\tilde{\phi}_{\alpha}^{+}(k)$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{V}}\int ds\;\exp(+ik\cdot
s)\;\tilde{\psi}_{\alpha}^{+}(s)$ (198)

where the $t$ or $\beta$ dependence is implicit.

##  11 Appendix B - Stochastic and QCF Equations

In this Appendix the details are set out of applying stochastic averaging to
give first order changes in the QCF describing the size of one Cooper pair and
the position correlations for two Cooper pairs.

###  11.1 One Cooper Pair - Time Evolution

For the case of time evolution the QCF at time $t+\delta t$ is related to that
at time $t$ via

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1},t+\delta
t)\,\widetilde{\phi}_{u}(k_{2},t+\delta
t)\;\widetilde{\phi}_{u}^{+}(k_{3},t+\delta
t)\,\widetilde{\phi}_{d}^{+}(k_{4},t+\delta t)}$ (203) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\mathop{\displaystyle\sum}\limits_{l_{1},q_{1}}(F_{d,d}(k_{1},l_{1};q_{1},\delta
t)\,\widetilde{\phi}_{d}(l_{1},t)+F_{d,u}(k_{1},l_{1};q_{1},\delta
t)\,\widetilde{\phi}_{u}(l_{1},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{2},q_{2}}(F_{u,d}(k_{2},l_{2};q_{2},\delta
t)\,\widetilde{\phi}_{d}(l_{2},t)+F_{u,u}(k_{2},l_{2};q_{2},\delta
t)\,\widetilde{\phi}_{u}(l_{2},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{3},q_{3}}(F_{u,d}^{+}(k_{3},l_{3};q_{3},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{3},t)+F_{u,u}^{+}(k_{3},l_{3};q_{3},\delta
t)\,\widetilde{\phi}_{u}^{+}(l_{3},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{4},q_{4}}(F_{d,d}^{+}(k_{4},l_{4};q_{4},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)+F_{d,u}^{+}(k_{4},l_{4};q_{4},\delta
t)\,\widetilde{\phi}_{u}^{+}(l_{4},t))\end{array}\right]_{StochAver}$

after substituting from Eq.(LABEL:Eq.ItoSMtmTimeEvln). Apart from the sums
over $l_{1,}q_{1},...,l_{4},q_{4}$ there are $16$ terms for which the
stochastic averages have to be evaluated. If the two terms in each of the
first $2$ factors are listed as $1,2$ and the two terms in each of the second
$2$ factors are listed as $a,b$ then the $16$ terms can be listed as
$11aa,12aa,21aa,22aa,11ab,12ab,21ab,22ab,....,22bb$.

In Ito stochastic equations we then use the factorisation result that depends
on the $F$, $F^{+}$ involving times later than $t$ (see Eq. (9.43) in Ref.
[19]

$\displaystyle\left[\begin{array}[]{c}F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta
t)\,\widetilde{\phi}_{\beta}(l_{a},t)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta
t)\,\widetilde{\phi}_{\delta}(l_{b},t)\\\
\times\,F_{\epsilon,\eta}^{+}(k_{c},l_{c};q_{c},\delta
t)\,\widetilde{\phi}_{\eta}^{+}(l_{c},t)\,F_{\theta,\kappa}^{+}(k_{d},l_{d};q_{d},\delta
t)\,\widetilde{\phi}_{\kappa}^{+}(l_{d},t)\end{array}\right]_{Stoch\,Aver}$
(208) $\displaystyle=$
$\displaystyle\overline{F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta
t)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta
t)\,F_{\epsilon,\eta}^{+}(k_{c},l_{c};q_{c},\delta
t)\,F_{\theta,\kappa}^{+}(k_{d},l_{d};q_{d},\delta t)}$
$\displaystyle\times\overline{\,\widetilde{\phi}_{\beta}(l_{a},t)\,\,\widetilde{\phi}_{\delta}(l_{b},t)\,\widetilde{\phi}_{\eta}^{+}(l_{c},t)\,\widetilde{\phi}_{\kappa}^{+}(l_{d},t)}$

We wish to evaluate the various
$\overline{F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta
t)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta
t)\,F_{\epsilon,\eta}^{+}(k_{c},l_{c};q_{c},\delta
t)\,F_{\theta,\kappa}^{+}(k_{d},l_{d};q_{d},\delta t)}$ correct to order
$\delta t$. We note that the $F_{d,d},F_{u,u},F_{d,d}^{+},F_{u,u}^{+}$ factors
are non-stochastic and linear in $\delta t$. Also the
$F_{d,u},F_{u,d},F_{d,u}^{+},F_{u,d}^{+}$ factors are stochastic and linear in
the stochastic Wiener increments
$\delta\widetilde{\omega}_{u,d}^{q},\delta\widetilde{\omega}_{d,u}^{q},\delta\widetilde{\omega}_{u+,d+}^{q},\delta\widetilde{\omega}_{d+,u+}^{q}$.
From Eq.(LABEL:Eq.PropsWienerIncrem) the stochastic average of a factor
involving an odd number $1,3$ of the
$\delta\widetilde{\omega}_{u,d}^{q},\delta\widetilde{\omega}_{d,u}^{q},\delta\widetilde{\omega}_{u+,d+}^{q},\delta\widetilde{\omega}_{d+,u+}^{q}$
with any non-stochastic factor is zero, if there are $2$ of these increments
the stochastic average is proportional to $\delta t$ and if there are $4$ of
these increments the stochastic average is proportional to $\delta
t^{2}$(apart from cases where there is zero Kronecka delta $\delta_{{\small
ab}}$ involved). For the case of two stochastic Wiener increments we note from
Eq. (LABEL:Eq.PropsWienerIncrem) that the non-zero stochastic averages of the
products are

$\displaystyle\overline{\delta\widetilde{\omega}_{u,d}^{q}\,\delta\widetilde{\omega}_{u,d}^{q}}$
$\displaystyle=$
$\displaystyle\overline{\delta\widetilde{\omega}_{d,u}^{q}\,\delta\widetilde{\omega}_{d,u}^{q}}=\delta
t$
$\displaystyle\overline{\delta\widetilde{\omega}_{u+,d+}^{q}\,\delta\widetilde{\omega}_{u+,d+}^{q}}$
$\displaystyle=$
$\displaystyle\overline{\delta\widetilde{\omega}_{d+,u+}^{q}\,\delta\widetilde{\omega}_{d+,u+}^{q}}=\delta
t$ (209)

With these considerations in mind, the non-zero cases are as follows:

$\displaystyle 12ba$ $\displaystyle:$
$\displaystyle\overline{F_{d,d}(k_{1},l_{1};q_{1},\delta
t)F_{u,u}(k_{2},l_{2};q_{2},\delta t)F_{u,u}^{+}(k_{3},l_{3};q_{3},\delta
t)\,F_{d,d}^{+}(k_{4},l_{4};q_{4},\delta t)\,}$ (210) $\displaystyle=$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}\right\\}\delta
t\right)\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{2}^{2}}{2m}\right\\}\delta
t\right)$
$\displaystyle\times\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{3}^{2}}{2m}\right\\}\delta
t\right)\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$ $\displaystyle\approx$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\,\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\,\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}$
$\displaystyle\times\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}-\frac{\hbar^{2}k_{3}^{2}}{2m}-\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$

where only terms linear in $\delta t$ have been retained.

$\displaystyle 12ab$ $\displaystyle:$
$\displaystyle\overline{F_{d,d}(k_{1},l_{1};q_{1},\delta
t)F_{u,u}(k_{2},l_{2};q_{2},\delta t)F_{u,d}^{+}(k_{3},l_{3};q_{3},\delta
t)\,F_{d,u}^{+}(k_{4},l_{4};q_{4},\delta t)\,}$ (211) $\displaystyle=$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}\right\\}\delta
t\right)\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{2}^{2}}{2m}\right\\}\delta
t\right)$
$\displaystyle\times\overline{\lambda\;\delta_{(k_{3}+q_{3}),l_{3}}\;\left\\{\delta\widetilde{\omega}_{u+,d+}^{q_{3}}+i\delta\widetilde{\omega}_{d+,u+}^{q_{3}}\right\\}\,\lambda\;\delta_{(k_{4}-q_{4}),l_{4}}\;\left\\{-\delta\widetilde{\omega}_{u+,d+}^{q_{4}}+i\delta\widetilde{\omega}_{d+,u+}^{q_{4}}\right\\}}$
$\displaystyle=$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\,\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\;\lambda\;\delta_{(k_{3}+q_{3}),l_{3}}\;\lambda\;\delta_{(k_{4}-q_{4}),l_{4}}$
$\displaystyle\times\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}\right\\}\delta
t\right)\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{2}^{2}}{2m}\right\\}\delta
t\right)$ $\displaystyle\times\left\\{-\delta_{q_{3},q_{4}}\delta
t-\delta_{q_{3},q_{4}}\delta t\right\\}$ $\displaystyle\approx$
$\displaystyle-2\lambda^{2}\,\\{\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\,\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\;\delta_{(k_{3}+q_{3}),l_{3}}\;\delta_{(k_{4}-q_{4}),l_{4}}\,\delta_{q_{3},q_{4}}\\}\,\delta
t$

where the Wiener increment properties in Eq.(LABEL:Eq.PropsWienerIncrem) has
been used and only terms linear in $\delta t$ have been retained.

$\displaystyle 21ba$ $\displaystyle:$
$\displaystyle\overline{F_{d,u}(k_{1},l_{1};q_{1},\delta
t)F_{u,d}(k_{2},l_{2};q_{2},\delta t)F_{u,u}^{+}(k_{3},l_{3};q_{3},\delta
t)\,F_{d,d}^{+}(k_{4},l_{4};q_{4},\delta t)\,}$ (212) $\displaystyle=$
$\displaystyle\overline{(+\lambda\;\delta_{(k_{1}+q_{1}),l_{1}}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q_{1}}-i\delta\widetilde{\omega}_{d,u}^{q_{1}}\right\\})\;(+\lambda\;\delta_{(k_{2}-q_{2}),l_{2}}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q_{2}}+i\delta\widetilde{\omega}_{d,u}^{q_{2}}\right\\}}$
$\displaystyle\times\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{3}^{2}}{2m}\right\\}\delta
t\right)\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$ $\displaystyle=$
$\displaystyle\lambda\;\delta_{(k_{1}+q_{1}),l_{1}}\,\lambda\;\delta_{(k_{2}-q_{2}),l_{2}}\;\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}$
$\displaystyle\times\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{3}^{2}}{2m}\right\\}\delta
t\right)\left(1+\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$ $\displaystyle\times\\{\delta_{q_{1},q_{2}}\delta
t+\delta_{q_{1},q_{2}}\delta t\\}$ $\displaystyle\approx$ $\displaystyle
2\lambda^{2}\,\\{\delta_{(k_{1}+q_{1}),l_{1}}\,\delta_{(k_{2}-q_{2}),l_{2}}\;\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}\,\delta_{q_{1},q_{2}}\\}\delta
t$

where only terms linear in $\delta t$ have been retained.

Substituting the results in Eqs.(210), (211) and (212) into Eq.
(LABEL:Eq.QCFTimeEvnStep1) and performing the sums over
$q_{1},q_{2},q_{3},q_{4}$ gives

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1},t+\delta
t)\,\widetilde{\phi}_{u}(k_{2},t+\delta
t)\;\widetilde{\phi}_{u}^{+}(k_{3},t+\delta
t)\,\widetilde{\phi}_{d}^{+}(k_{4},t+\delta t)}$ $\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}-\frac{\hbar^{2}k_{3}^{2}}{2m}-\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$
$\displaystyle-2\lambda^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}\\}\,\delta
t\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{d}^{+}(l_{3},t)\,\widetilde{\phi}_{u}^{+}(l_{4},t)}$
$\displaystyle+2\lambda^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{(k_{1}+k_{2}),(l_{1}+l_{2})}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\\}\,\delta
t\times\overline{\widetilde{\phi}_{u}(l_{1},t)\,\widetilde{\phi}_{d}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\left(1-\frac{i}{\hbar}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}-\frac{\hbar^{2}k_{3}^{2}}{2m}-\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta
t\right)$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$
$\displaystyle+2\lambda^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}\\}\,\delta
t\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$
$\displaystyle-2\lambda^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{(k_{1}+k_{2}),(l_{1}+l_{2})}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\\}\,\delta
t\times\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{u}^{+}(l_{3},t)\,\widetilde{\phi}_{d}^{+}(l_{4},t)}$

where we have used
$\mathop{\displaystyle\sum}\limits_{q_{3}q_{4}}\delta_{(k_{3}+q_{3}),l_{3}}\;\delta_{(k_{4}-q_{4}),l_{4}}\,\delta_{q_{3},q_{4}}=\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}$
and $\mathop{\displaystyle\sum}\limits_{q_{1}}\delta_{q_{1},0}=1$. Also, in
the second last line we have interchanged $l_{3}$ and $l_{4}$ and used the
anti-commuting feature of the Grassmann stochastic momentum fields
$\widetilde{\phi}_{d}^{+}(l_{3},t)$ and $\,\widetilde{\phi}_{u}^{+}(l_{4},t)$
to place these in the opposite order in the stochastic average
$\overline{\widetilde{\phi}_{d}(l_{1},t)\,\widetilde{\phi}_{u}(l_{2},t)\;\widetilde{\phi}_{d}^{+}(l_{3},t)\,\widetilde{\phi}_{u}^{+}(l_{4},t)}$,
with similar steps for the last line.

###  11.2 One Cooper Pair - Temperature Evolution

For the case of temperature evolution the QCF at time $\beta+\delta\beta$ is
related to that at time $t$ via

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1},\beta+\delta\beta)\,\widetilde{\phi}_{u}(k_{2},\beta+\delta\beta)\;\widetilde{\phi}_{u}^{+}(k_{3},\beta+\delta\beta)\,\widetilde{\phi}_{d}^{+}(k_{4},\beta+\delta\beta)}$
(218) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\mathop{\displaystyle\sum}\limits_{l_{1},q_{1}}(F_{d,d}(k_{1},l_{1};q_{1},\delta\beta)\,\widetilde{\phi}_{d}(l_{1},\beta)+F_{d,u}(k_{1},l_{1};q_{1},\delta\beta)\,\widetilde{\phi}_{u}(l_{1},\beta))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{2},q_{2}}(F_{u,d}(k_{2},l_{2};q_{2},\delta\beta)\,\widetilde{\phi}_{d}(l_{2},\beta)+F_{u,u}(k_{2},l_{2};q_{2},\delta\beta)\,\widetilde{\phi}_{u}(l_{2},\beta))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{3},q_{3}}(F_{u,d}^{+}(k_{3},l_{3};q_{3},\delta\beta)\,\widetilde{\phi}_{d}^{+}(l_{3},\beta)+F_{u,u}^{+}(k_{3},l_{3};q_{3},\delta\beta)\,\widetilde{\phi}_{u}^{+}(l_{3},\beta))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{4},q_{4}}(F_{d,d}^{+}(k_{4},l_{4};q_{4},\delta\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)+F_{d,u}^{+}(k_{4},l_{4};q_{4},\delta\beta)\,\widetilde{\phi}_{u}^{+}(l_{4},\beta))\end{array}\right]_{StochAver}$

after substituting from Eq.(LABEL:Eq.ItoSMtmFTemper). Apart from the sums over
$l_{1,}q_{1},...,l_{4},q_{4}$ there are $16$ terms for which the stochastic
averages have to be evaluated. .If the two terms in each of the first $2$
factors are listed as $1,2$ and the two terms in each of the second $2$
factors are listed as $a,b$ then the $16$ terms can be listed as
$11aa,12aa,21aa,22aa,11ab,12ab,21ab,22ab,....,22bb$.

In Ito stochastic equations we then use the factorisation result that depends
on the $F$, $F^{+}$ involving times later than $t$ (see Eq. (9.43) in Ref.
[19]

$\displaystyle\left[\begin{array}[]{c}F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta\beta)\,\widetilde{\phi}_{\beta}(l_{1},\beta)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta\beta)\,\widetilde{\phi}_{\delta}(l_{b},\beta)\\\
\times\,F_{\epsilon,\eta}^{+}(k_{c},l_{c};q_{c},\delta\beta)\,\widetilde{\phi}_{\eta}^{+}(l_{c},\beta)\,F_{\theta,\kappa}^{+}(k_{d},l_{d};q_{d},\delta\beta)\,\widetilde{\phi}_{\kappa}^{+}(l_{d},\beta)\end{array}\right]_{Stoch\,Aver}$
(223) $\displaystyle=$
$\displaystyle\overline{F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta\beta)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta\beta)\,F_{\epsilon,\eta}^{+}(k_{c},l_{c};q_{c},\delta\beta)\,F_{\theta,\kappa}^{+}(k_{d},l_{d};q_{d},\delta\beta)}$
$\displaystyle\times\overline{\,\widetilde{\phi}_{\beta}(l_{a},\beta)\,\,\widetilde{\phi}_{\delta}(l_{b},\beta)\,\widetilde{\phi}_{\eta}^{+}(l_{c},\beta)\,\widetilde{\phi}_{\kappa}^{+}(l_{d},\beta)}$

We wish to evaluate the various
$\overline{F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta\beta)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta\beta)\,F_{\epsilon,\eta}^{+}(k_{c},l_{c};q_{c},\delta\beta)\,F_{\theta,\kappa}^{+}(k_{d},l_{d};q_{d},\delta\beta)}$
correct to order $\delta\beta$. We note that the
$F_{d,d},F_{u,u},F_{d,d}^{+},F_{u,u}^{+}$ factors are non-stochastic and
linear in $\delta\beta$. Also the $F_{d,u},F_{u,d},F_{d,u}^{+},F_{u,d}^{+}$
factors are stochastic and linear in the stochastic Wiener increments
$\delta\widetilde{\omega}_{u,d}^{q},\delta\widetilde{\omega}_{d,u}^{q},\delta\widetilde{\omega}_{u+,d+}^{q},\delta\widetilde{\omega}_{d+,u+}^{q}$.
From Eq.(LABEL:Eq.PropsWienerIncrem) the stochastic average of a factor
involving an odd number $1,3$ of the
$\delta\widetilde{\omega}_{u,d}^{q},\delta\widetilde{\omega}_{d,u}^{q},\delta\widetilde{\omega}_{u+,d+}^{q},\delta\widetilde{\omega}_{d+,u+}^{q}$
with any non-stochastic factor is zero, if there are $2$ of these increments
the stochastic average is proportional to $\delta\beta$ and if there are $4$
of these increments the stochastic average is proportional to
$\delta\beta^{2}$(apart from cases where there is zero Kronecka delta
$\delta_{{\small ab}}$ involved). For the case of two stochastic Wiener
increments we note from Eq. (LABEL:Eq.PropsWienerIncrem) that the non-zero
stochastic averages of the products are

$\displaystyle\overline{\delta\widetilde{\omega}_{u,d}^{q}\,\delta\widetilde{\omega}_{u,d}^{q}}$
$\displaystyle=$
$\displaystyle\overline{\delta\widetilde{\omega}_{d,u}^{q}\,\delta\widetilde{\omega}_{d,u}^{q}}=\delta\beta$
$\displaystyle\overline{\delta\widetilde{\omega}_{u+,d+}^{q}\,\delta\widetilde{\omega}_{u+,d+}^{q}}$
$\displaystyle=$
$\displaystyle\overline{\delta\widetilde{\omega}_{d+,u+}^{q}\,\delta\widetilde{\omega}_{d+,u+}^{q}}=\delta\beta$
(224)

With these considerations in mind, the non-zero cases are as follows:

$\displaystyle 12ba$ $\displaystyle:$
$\displaystyle\overline{F_{d,d}(k_{1},l_{1};q_{1},\delta\beta)F_{u,u}(k_{2},l_{2};q_{2},\delta\beta)F_{u,u}^{+}(k_{3},l_{3};q_{3},\delta\beta)\,F_{d,d}^{+}(k_{4},l_{4};q_{4},\delta\beta)\,}$
(225) $\displaystyle=$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}\right\\}\delta\beta\right)\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{2}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{3}^{2}}{2m}\right\\}\delta\beta\right)\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\approx$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\,\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\,\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}$
$\displaystyle\times\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}+\frac{\hbar^{2}k_{3}^{2}}{2m}+\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$

where only terms linear in $\delta t$ have been retained.

$\displaystyle 12ab$ $\displaystyle:$
$\displaystyle\overline{F_{d,d}(k_{1},l_{1};q_{1},\delta\beta)F_{u,u}(k_{2},l_{2};q_{2},\delta\beta)F_{u,d}^{+}(k_{3},l_{3};q_{3},\delta\beta)\,F_{d,u}^{+}(k_{4},l_{4};q_{4},\delta\beta)\,}$
(226) $\displaystyle=$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}\right\\}\delta\beta\right)\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{2}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\overline{\eta\;\delta_{(k_{3}+q_{3}),l_{3}}\;\left\\{\delta\widetilde{\omega}_{u+,d+}^{q_{3}}+i\delta\widetilde{\omega}_{d+,u+}^{q_{3}}\right\\}\,\eta\;\delta_{(k_{4}-q_{4}),l_{4}}\;\left\\{\delta\widetilde{\omega}_{u+,d+}^{q_{4}}-i\delta\widetilde{\omega}_{d+,u+}^{q_{4}}\right\\}}$
$\displaystyle=$
$\displaystyle\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\,\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\;\eta\;\delta_{(k_{3}+q_{3}),l_{3}}\;\eta\;\delta_{(k_{4}-q_{4}),l_{4}}$
$\displaystyle\times\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}\right\\}\delta\beta\right)\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{2}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\left\\{+\delta_{q_{3},q_{4}}\delta\beta+\delta_{q_{3},q_{4}}\delta\beta\right\\}$
$\displaystyle\approx$
$\displaystyle+2\eta^{2}\,\\{\delta_{q_{1},0}\;\delta_{k_{1},l_{1}}\,\delta_{q_{2},0}\;\delta_{k_{2},l_{2}}\;\delta_{(k_{3}+q_{3}),l_{3}}\;\delta_{(k_{4}-q_{4}),l_{4}}\,\delta_{q_{3},q_{4}}\\}\,\delta\beta$

where the Wiener increment properties in Eq.(LABEL:Eq.PropsWienerIncrem) has
been used and only terms linear in $\delta t$ have been retained.

$\displaystyle 21ba$ $\displaystyle:$
$\displaystyle\overline{F_{d,u}(k_{1},l_{1};q_{1},\delta\beta)F_{u,d}(k_{2},l_{2};q_{2},\delta\beta)F_{u,u}^{+}(k_{3},l_{3};q_{3},\delta\beta)\,F_{d,d}^{+}(k_{4},l_{4};q_{4},\delta\beta)\,}$
(227) $\displaystyle=$
$\displaystyle\overline{(+\eta\;\delta_{(k_{1}+q_{1}),l_{1}}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q_{1}}-i\delta\widetilde{\omega}_{d,u}^{q_{1}}\right\\})\;(+\eta\;\delta_{(k_{2}-q_{2}),l_{2}}\;\left\\{\delta\widetilde{\omega}_{u,d}^{q_{2}}+i\delta\widetilde{\omega}_{d,u}^{q_{2}}\right\\}}$
$\displaystyle\times\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{3}^{2}}{2m}\right\\}\delta\beta\right)\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle=$
$\displaystyle\eta\;\delta_{(k_{1}+q_{1}),l_{1}}\,\eta\;\delta_{(k_{2}-q_{2}),l_{2}}\;\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}$
$\displaystyle\times\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{3}^{2}}{2m}\right\\}\delta\beta\right)\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\\{\delta_{q_{1},q_{2}}\delta\beta+\delta_{q_{1},q_{2}}\delta\beta\\}$
$\displaystyle\approx$ $\displaystyle
2\eta^{2}\,\\{\delta_{(k_{1}+q_{1}),l_{1}}\,\delta_{(k_{2}-q_{2}),l_{2}}\;\delta_{q_{3},0}\;\delta_{k_{3},l_{3}}\,\delta_{q_{4},0}\;\delta_{k_{4},l_{4}}\,\delta_{q_{1},q_{2}}\\}\delta\beta$

where only terms linear in $\delta\beta$ have been retained.

Substituting the results in Eqs.(225), (226) and (227) into Eq.
(LABEL:Eq.QCFTempEvnStep1) and performing the sums over
$q_{1},q_{2},q_{3},q_{4}$ gives

$\displaystyle\overline{\widetilde{\phi}_{d}(k_{1},\beta+\delta\beta)\,\widetilde{\phi}_{u}(k_{2},\beta+\delta\beta)\;\widetilde{\phi}_{u}^{+}(k_{3},\beta+\delta\beta)\,\widetilde{\phi}_{d}^{+}(k_{4},\beta+\delta\beta)}$
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}+\frac{\hbar^{2}k_{3}^{2}}{2m}+\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$
$\displaystyle+2\eta^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}\\}\,\delta\beta\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{d}^{+}(l_{3},\beta)\,\widetilde{\phi}_{u}^{+}(l_{4},\beta)}$
$\displaystyle+2\eta^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{(k_{1}+k_{2}),(l_{1}+l_{2})}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\\}\,\delta\beta\times\overline{\widetilde{\phi}_{u}(l_{1},\beta)\,\widetilde{\phi}_{d}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$
$\displaystyle=$
$\displaystyle\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\left(1-\frac{1}{2}\left\\{\frac{\hbar^{2}k_{1}^{2}}{2m}+\frac{\hbar^{2}k_{2}^{2}}{2m}+\frac{\hbar^{2}k_{3}^{2}}{2m}+\frac{\hbar^{2}k_{4}^{2}}{2m}\right\\}\delta\beta\right)$
$\displaystyle\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$
$\displaystyle-2\eta^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{k_{1},l_{1}}\,\delta_{k_{2},l_{2}}\,\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}\\}\,\delta\beta\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$
$\displaystyle-2\eta^{2}\,\mathop{\displaystyle\sum}\limits_{l_{1}l_{2}l_{3}l_{4}}\\{\delta_{(k_{1}+k_{2}),(l_{1}+l_{2})}\,\delta_{k_{3},l_{3}}\,\delta_{k_{4},l_{4}}\\}\,\delta\beta\times\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{u}^{+}(l_{3},\beta)\,\widetilde{\phi}_{d}^{+}(l_{4},\beta)}$

where we have used
$\mathop{\displaystyle\sum}\limits_{q_{3}q_{4}}\delta_{(k_{3}+q_{3}),l_{3}}\;\delta_{(k_{4}-q_{4}),l_{4}}\,\delta_{q_{3},q_{4}}=\delta_{(k_{3}+k_{4}),(l_{3}+l_{4})}$
and $\mathop{\displaystyle\sum}\limits_{q_{1}}\delta_{q_{1},0}=1$. Also, in
the second last line we have interchanged $l_{3}$ and $l_{4}$ and used the
anti-commuting feature of the Grassmann stochastic momentum fields
$\widetilde{\phi}_{d}^{+}(l_{3},\beta)$ and
$\,\widetilde{\phi}_{u}^{+}(l_{4},\beta)$ to place these in the opposite order
in the stochastic average
$\overline{\widetilde{\phi}_{d}(l_{1},\beta)\,\widetilde{\phi}_{u}(l_{2},\beta)\;\widetilde{\phi}_{d}^{+}(l_{3},\beta)\,\widetilde{\phi}_{u}^{+}(l_{4},\beta)}$,
with similar steps for the last line.

###  11.3 Two Cooper Pairs - Time Evolution

For the case of time evolution the QCF for two Cooper pairs at time $t+\delta
t$ is related to that at time $t$ via

$\displaystyle\left[\begin{array}[]{c}\widetilde{\phi}_{d}(k_{1},t+\delta
t)\,\widetilde{\phi}_{d}(k_{2},t+\delta
t)\;\widetilde{\phi}_{u}(k_{3},t+\delta
t)\,\widetilde{\phi}_{u}(k_{4},t+\delta t)\\\
\times\;\widetilde{\phi}_{u}^{+}(k_{5},t+\delta
t)\,\widetilde{\phi}_{u}^{+}(k_{6},t+\delta
t)\widetilde{\phi}_{d}^{+}(k_{7},t+\delta
t)\,\widetilde{\phi}_{d}^{+}(k_{8},t+\delta t)\end{array}\right]_{StochAver}$
(231) $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\mathop{\displaystyle\sum}\limits_{l_{1},q_{1}}(F_{d,d}(k_{1},l_{1};q_{1},\delta
t)\,\widetilde{\phi}_{d}(l_{1},t)+F_{d,u}(k_{1},l_{1};q_{1},\delta
t)\,\widetilde{\phi}_{u}(l_{1},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{2},q_{2}}(F_{d,d}(k_{2},l_{2};q_{2},\delta
t)\,\widetilde{\phi}_{d}(l_{2},t)+F_{d,u}(k_{2},l_{2};q_{2},\delta
t)\,\widetilde{\phi}_{u}(l_{2},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{3},q_{3}}(F_{u,d}(k_{3},l_{3};q_{3},\delta
t)\,\widetilde{\phi}_{d}(l_{3},t)+F_{u,u}(k_{3},l_{3};q_{3},\delta
t)\,\widetilde{\phi}_{u}(l_{3},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{4},q_{4}}(F_{u,d}(k_{4},l_{4};q_{4},\delta
t)\,\widetilde{\phi}_{d}(l_{4},t)+F_{u,u}(k_{4},l_{4};q_{4},\delta
t)\,\widetilde{\phi}_{u}(l_{4},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{5},q_{5}}(F_{u,d}^{+}(k_{5},l_{5};q_{5},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{5},t)+F_{u,u}^{+}(k_{5},l_{5};q_{5},\delta
t)\,\widetilde{\phi}_{u}^{+}(l_{5},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{6},q_{6}}(F_{u,d}^{+}(k_{6},l_{6};q_{6},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{6},t)+F_{u,u}^{+}(k_{6},l_{6};q_{6},\delta
t)\,\widetilde{\phi}_{u}^{+}(l_{6},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{7},q_{7}}(F_{d,d}^{+}(k_{7},l_{7};q_{7},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{7},t)+F_{d,u}^{+}(k_{7},l_{7};q_{7},\delta
t)\,\widetilde{\phi}_{u}^{+}(l_{7},t))\\\
\times\mathop{\displaystyle\sum}\limits_{l_{8},q_{8}}(F_{d,d}^{+}(k_{8},l_{8};q_{8},\delta
t)\,\widetilde{\phi}_{d}^{+}(l_{8},t)+F_{d,u}^{+}(k_{8},l_{8};q_{8},\delta
t)\,\widetilde{\phi}_{u}^{+}(l_{8},t))\end{array}\right]_{StochAver}$ (240)

after substituting from Eq.(LABEL:Eq.ItoSMtmTimeEvln). Apart from the sums
over $l_{1,}q_{1},...,l_{8},q_{8}$ there are $256$ terms for which the
stochastic averages have to be evaluated. If the two terms in each of the
first $4$ factors are listed as $1,2$ and the two terms in each of the second
$4$ factors are listed as $a,b$ then the $256$ terms can be listed as

$1111aaaa,1211aaaa,2111aaaa,2211aaaa,....,2222aaaa,$

$1111abaa,1211abaa,2111abaa,2211abaa,....,2222abaa,....,1111bbbb,....2222bbbb.$.

In Ito stochastic equations we then use the factorisation result that depends
on the $F$, $F^{+}$ involving times later than $t$ (see Eq. (9.43) in Ref.
[19]

$\displaystyle\left[\begin{array}[]{c}F_{\alpha,\beta}(k_{a},l_{a};q_{a},\delta
t)\,\widetilde{\phi}_{\beta}(l_{a},t)\,F_{\gamma,\delta}(k_{b},l_{b};q_{b},\delta
t)\,\widetilde{\phi}_{\delta}(l_{b},t)\\\ \times
F_{\epsilon,\eta}(k_{c},l_{c};q_{c},\delta
t)\,\widetilde{\phi}_{\eta}(l_{c},t)\,F_{\theta,\kappa}(k_{d},l_{d};q_{d},\delta
t)\,\widetilde{\phi}_{\kappa}(l_{d},t)\\\ \times
F_{\pi,\xi}^{+}(k_{e},l_{e};q_{e},\delta
t)\,\widetilde{\phi}_{\xi}^{+}(l_{e},t)\,F_{\zeta,\tau}^{+}(k_{f},l_{f};q_{f},\delta
t)\,\widetilde{\phi}_{\tau}^{+}(l_{f},t)\\\