text:
# Accuracy and convergence of the backward Monte-Carlo method

## 1 Introduction

The simplest parabolic equation is the canonical diffusion equation,

$$\frac{f}{t}=\frac{}{x}D\frac{f}{x},x,0tT,$$
(1)
with some initial condition $`f(x,0)=\mathrm{\Phi }(x)`$. The conventional, forward, weighted Monte-Carlo solution to Eq. (1) is given by

$$f_{}(x,T)=N^1\underset{i=1}{\overset{N}{}}\mathrm{\Phi }(X_i^{}(0))\times \delta (xX_i^{}(T)),$$
(2)
where the stochastic variables $`X_i^{}(T)`$ are found by following the stochastic trajectories given by the forward Monte-Carlo difference equation of motion:

$$X_i^{}(t+\mathrm{\Delta }t)=X_i^{}(t)+\mu \mathrm{\Delta }t+\zeta \sigma \sqrt{\mathrm{\Delta }t},i=1,\mathrm{},N,$$
(3)
where $`\mu =D/x`$, $`\sigma =\sqrt{2D}`$, and $`\zeta `$ is a zero-mean, unit-variance Gaussian random number, $`\zeta N(0,1)`$. The points where the particles are launched, $`X_i^{}(0)`$, must be chosen so that $`f_{}(x,T=0)`$, as given by Eq. (2), approximates $`\mathrm{\Phi }(x)`$. The algorithm is illustrated in Fig. 1 below.

The recently introduced *backward* Monte-Carlo method  is both strikingly similar and fundamentally different from its older, forward sibling. The backward solution to Eq. (1) is

$$f_{}(x,T)=N^1\underset{i=1}{\overset{N}{}}\mathrm{\Phi }(X_i^{}(0)),$$
(4)
where the stochastic variables $`X_i^{}(0)`$ are found by following the stochastic trajectories given by the *backward* Monte-Carlo difference equation of motion:

$$X_i^{}(t\mathrm{\Delta }t)=X_i^{}(t)+\mu \mathrm{\Delta }t+\zeta \sigma \sqrt{\mathrm{\Delta }t},X_i^{}(T)=x,i=1,\mathrm{},N.$$
(5)
Comparing the forward and backward solutions, Eqs. (2) and (4), respectively, the appearance of $`\delta `$-functions in the forward solution is the most striking difference. To understand why the backward method produces a smooth solution, i.e. without $`\delta `$-functions, we have to go back to its derivation , which relies heavily on the Feynman-Kac formula . In fact, Eq. (4), supplemented by Eq. (5), is just the numerical approximation of the stochastic Feynman-Kac representation of Eq. (1). The Feynman-Kac formula relates a parabolic equation to its *naturally associated stochastic differential equation* (SDE). This SDE governs microscopic motion going *backward* in time, and Eq. (5) is its difference approximation. As a consequence of time running in different directions in the macroscopic and microscopic world, we get the natural initial condition, $`X_i^{}(T)=x`$ in Eq. (5); i.e. we launch the particles from the point in extended phase space where we want to know the solution to the parabolic equation and let them work their way backward in time to sample from the initial condition $`\mathrm{\Phi }(X_i^{}(0))`$ \[see Fig. 2 above\]. As is shown in appendix A, the forward Monte-Carlo method can be derived by forcing the naturally associated SDE to govern microscopic motion going *forward* in time. This coercion is responsible for both the $`\delta `$-functions in the forward Monte-Carlo solution, Eq. (2), and the “fuzzy initial condition”, $`X_i^{}(0)\{X_i^{}(0)|f_{}(x,T=0)\mathrm{\Phi }(x)\}`$, imposed on the forward Monte-Carlo equation of motion, Eq. (3) \[see appendix A for details\].

However, the backward algorithm was not developed for esthetic reasons. The motivation originally came from the frustration over the extremely inefficient use of test particles of the forward method in cases where only the solution in a small part of phase space contributes significantly to the physics of interest. An example of this type of problem is the heating of fusion plasmas by cyclotron-resonance absorption on ions, where only the high-energy tail of the solution is of relevance; but most test particles are wasted on the almost Maxwellian bulk, which has been simulated with e.g. the FIDO forward Monte-Carlo code . Earlier claims  that the backward method is vastly superior in these cases will be substantiated in section 3.

Before that, the sources of numerical error, for both the forward and the backward method, will be identified and briefly discussed in the next section below. In section 3 we will compare the forward and backward Monte-Carlo solutions to a simple model equation. In section 4, it will be shown how a minor modification of the backward Monte-Carlo difference equation of motion, Eq. (5), leads to improved time convergence of the backward solution, Eq. (4). Finally, a summary of the main findings of this article follows in section 5.

## 2 Sources of error

The greatest disadvantage of the Monte-Carlo method is the unavoidable statistical error caused by the use of random numbers. With the conventional, forward Monte-Carlo method, the $`\delta `$-functions in the solution \[see Eq. (2)\] make binning of the particles necessary. With the number of particles in each bin proportional to $`N/N_{\mathrm{bin}}`$, the statistical error becomes $`𝒪((N/N_{\mathrm{bin}})^{1/2})`$ with obvious notations. Backing out the solution by linear interpolation between the bins introduces an error $`𝒪(N_{\mathrm{bin}}^2)`$.

As was mentioned previously, the backward method was developed for cases where the solution is of interest only in a point or small part of phase space. In such situations it reduces the statistical error to $`𝒪(N^{1/2})`$. Even in cases where the backward method is used to calculate global solutions, and the statistical error becomes $`𝒪((N/N_{\mathrm{bin}})^{1/2})`$, it does not have a finite-bin-size error because the absence of $`\delta `$-functions in the solution \[see Eq. (4)\] makes binning superfluous.

The finite time step, $`\mathrm{\Delta }t`$, also introduces an error term. An analysis of this time-step error requires some care; some mathematical subtleties are involved. There are two ways through which the time-step error can enter the numerical solution.

First, directly through the numerical approximation of the Feynman-Kac formula itself. I.e, even if we knew the *exact* stochastic variables $`X_i(0)`$, the backward solution, Eq. (4), would have an error due to the finite time step. To estimate the magnitude of this error, we need to go back to the derivation of Itô’s formula in Appendix A of Ref.  where terms of order $`𝒪(\mathrm{\Delta }t^{3/2})`$ and higher were neglected. However, in the same way that the $`𝒪(\mathrm{\Delta }t^{1/2})`$-term, which was kept, becomes zero after averaging, so does the $`𝒪(\mathrm{\Delta }t^{3/2})`$-term. So, even if the $`X_i^{}(0)`$ were known to all orders, $`X_i^{}(0)X_i(0)`$, the backward solution Eq. (4), would have an error $`𝒪(\mathrm{\Delta }t^2)`$. The error of the forward solution, Eq. (4) would of course be of the same order.

The other way in which the finite time step introduces an error into the solution is through an error in $`X_i^{}(0)`$. Again, going back to Ref. , we see that the drift term, $`\mu \mathrm{\Delta }t`$, on the right-hand side of the backward Monte-Carlo equation of motion, Eq. (5), has an error $`𝒪(\mathrm{\Delta }t^{3/2})`$, and the diffusive term, $`\zeta \sigma \sqrt{\mathrm{\Delta }t}`$, has an error $`𝒪(\mathrm{\Delta }t)`$. The different character of these two error terms should be noted; the $`𝒪(\mathrm{\Delta }t^{3/2})`$-term is deterministic whereas the $`𝒪(\mathrm{\Delta }t)`$-term is stochastic. It can be shown (see e.g. Ref. ) that after averaging over all the $`X_i^{}(0)`$, the approximation Eq. (5) introduces an $`𝒪(\mathrm{\Delta }t)`$ error into the backward solution. Again, the situation is equivalent for the forward method.

In section 4 we will discuss how the time-step error of the backward solution can be reduced. But before that, we will study the backward and forward solutions to a simple model equation in the next section.

## 3 Validation and benchmarking

To validate and benchmark the backward method and to compare it to the conventional forward method, we will solve a simple model equation that has an analytic solution. We have chosen the Lorentz equation that has been used to model pitch-angle scattering in plasmas :

$$\frac{f}{t}=\frac{}{x}D\frac{f}{x},1x1,0tT,$$
(6)
where the diffusion coefficient $`D=(1+x)(1x)`$ and the initial condition is $`f(x,0)=_{\mathrm{}=0}^{\mathrm{}}(\mathrm{}+1/2)P_{\mathrm{}}(x_0)P_{\mathrm{}}(x)e^{\mathrm{}(\mathrm{}+1)T_0}`$, where $`P_{\mathrm{}}(x)`$ is the Legendre polynomial. The analytic solution is $`f(x,T)=_{\mathrm{}=0}^{\mathrm{}}(\mathrm{}+1/2)P_{\mathrm{}}(x_0)P_{\mathrm{}}(x)e^{\mathrm{}(\mathrm{}+1)(T_0+T)}`$.

In the following we will use the parameters: $`x_0=0.9`$, $`T_0=0.1`$, and $`T=0.1`$; while $`N`$, $`N_{\mathrm{bin}}`$, and $`\mathrm{\Delta }t`$ will be varied so that their impact on the error can be studied. Before we get into a detailed analysis of accuracy and convergence, we show the solutions (with $`N=2\times 10^5`$, $`N_{\mathrm{bin}}=20`$, and $`\mathrm{\Delta }t=10^2`$) to the Lorentz equation (6) in Fig. 3 below.

Solving for $`f(x,T)`$ over the whole interval $`1x1`$, both the backward, $`f_{}(x,T)`$, and the forward, $`f_{}(x,T)`$ Monte-Carlo solutions are indistinguishable from the exact solution. However, solving only in the subinterval $`0.98x1`$, the backward solution is a much better approximation of the exact solution due to its much smaller statistical error. Note especially that the forward solution drops to zero for $`x0.993`$ because not a single particle finds its way into the last few bins \[compare also with the sketch in Fig. 1\]. The ability to focus in on a small subinterval and calculate efficiently (i.e. without wasting the vast majority of the particles) the local solution is one of the main strengths of the backward Monte-Carlo method.

For a detailed error analysis we need a well-defined measure of the error. We will use the local relative error $`\epsilon {}_{}{}^{}=|f_{}(x_0,T)f(x_0,T)|/f(x_0,T)`$. In the following we will plot the logarithm of the error against the logarithm of one of the parameters $`N`$, $`N_{\mathrm{bin}}`$, and $`\mathrm{\Delta }t`$, while keeping the other two parameters fixed. First we compare the statistical error of the forward and the backward solutions. To single out the statistical error from the finite bin-size and time-step errors we make the latter two small by choosing $`N_{\mathrm{bin}}=20`$, and $`\mathrm{\Delta }t=10^4`$. With $`N`$ spanning six decades ($`N=1,2,5,10,20,50,\mathrm{},10^6`$), we plot the error in the left frame of Fig. 4 below.

The dotted lines are least-square fits whose slopes approximate the exponents that determine the scaling of the error. With the forward method, the bins adjacent to $`x_0`$ are empty for small $`N`$ ($`N=1,2,20`$). As a result, $`\epsilon {}_{}{}^{}=1`$ for these values of $`N`$, and the slope becomes flat. The low-$`N`$ data points have thus been excluded from the least-square fit to $`\mathrm{log}\epsilon _{}`$. As can be seen, the statistical error \[$`𝒪(N^{1/2})`$ with $`N_{\mathrm{bin}}`$ fixed\] then scales as predicted. The exponents -0.43 and -0.41 are as close to the theoretical value of -1/2 as could be expected, given the fact that the other small but non-zero error terms always tend to flatten the slope. The forward and backward methods thus converge at the same rate as $`N`$ is increased, but the forward method needs $`N_{\mathrm{bin}}`$ times more particles to achieve the same accuracy as the backward method. This can easily translate into orders of magnitude in terms of execution time.

Next, we study the bin-size error (and the statistical error) by varying $`N_{\mathrm{bin}}`$ over four decades ($`N_{\mathrm{bin}}=1,2,5,\mathrm{},10^4`$) while keeping fixed $`N=10^4`$ and $`\mathrm{\Delta }t=10^4`$. The result is shown in the right frame of Fig. 4 above. The first observation is that the backward solution is completely unaffected by the bin size; $`\epsilon _{}`$ is constant and the slope is zero. The error in the forward solution, $`\epsilon _{}`$, exhibits a more interesting behavior with different scalings for small and large $`N_{\mathrm{bin}}`$. With few bins there are many particles in each bin, so the statistical error is small. The bin-size error, $`𝒪(N_{\mathrm{bin}}^2)`$, caused by backing out the solution by linear interpolation between these huge bins, however, becomes large. For small $`N_{\mathrm{bin}}`$, the slope of $`\mathrm{log}\epsilon _{}`$ is -1.90, in excellent agreement with the theoretical value of -2. With many bins, the bin-size error becomes small, but with few particles in each bin, the statistical error \[$`𝒪(N_{\mathrm{bin}}^{1/2})`$ with $`N`$ fixed\] becomes dominant as is evidenced by the slope ($`0.481/2`$) for large $`N_{\mathrm{bin}}`$. The shape of the $`\mathrm{log}\epsilon _{}`$ graph illustrates the conflict between resolution and statistics that is inherent for the forward Monte-Carlo method. The backward solution, however, can be calculated in two points arbitrarily close without affecting the statistical error. This is yet another of the main strengths of the backward method.

Turning now to the time-step error, we fix $`N=2\times 10^5`$ and $`N_{\mathrm{bin}}=20`$ and vary $`\mathrm{\Delta }t=10^5,2\times 10^5,5\times 10^5,\mathrm{},10^1`$. As can be seen in Fig. 5 below,

both the forward and the backward methods converge as predicted (the respective slopes are $`1`$) down to a time step $`\mathrm{\Delta }t10^3`$. At this point the statistical error becomes dominant.

## 4 A higher-order backward method

As was discussed in section 2, the dominant $`𝒪(\mathrm{\Delta }t)`$ time-step error in the Monte-Carlo equations of motion, Eqs. (3) and (5), comes from the diffusive term $`\zeta \sigma \sqrt{\mathrm{\Delta }t}`$. There is thus reason to suspect that a higher-order diffusive term might result in better overall time convergence, by reducing the time-step error to $`𝒪(\mathrm{\Delta }t^{3/2})`$. In appendix B it is shown that the next-higher-order Monte-Carlo equation of motion is

$$X_i^{}(t\pm \mathrm{\Delta }t)=X_i^{}(t)+\frac{1}{2}(1+\zeta ^2)\mu \mathrm{\Delta }t+\zeta \sigma \sqrt{\mathrm{\Delta }t}+𝒪(\mathrm{\Delta }t^{3/2}).$$
(7)
In the left frame of Fig. 6 below the time-convergence study presented in Fig. 5 above is repeated using the higher-order approximation of Eq. (7). In the right frame exactly the same experiment is repeated with ten times more particles to further reduce the statistical error.

Somewhat surprisingly, the forward solution does not converge at all. The backward solution, however, exhibits the faster time convergence we had hoped for (the slope is $`3/2`$).

Milshtein, who was not explicitly concerned with solving parabolic equations, did investigate the scaling of the error of expectation values similar to the ones of interest here . His results are consistent with the $`𝒪(\mathrm{\Delta }t)`$ error of the forward and backward lower-order solutions. However, for the higher-order methods introduced in this section, Milshtein’s results would still indicate a $`𝒪(\mathrm{\Delta }t)`$ error for both the forward and backward solutions, in clear disagreement with the scalings of Fig. 6. We hope to be able to resolve the discrepancies in future work.

## 5 Summary

We have shown that the backward Monte-Carlo method works as expected, i.e. it dramatically reduces the statistical error in situations where the solution is sought only in a small part of phase space. Furthermore, even in cases where we solve for the global solution, the backward method removes the conflict between resolution and statistics. This has great practical significance e.g. when the gradient of the solution is of interest.

We have also shown how a remarkably simple modification of the backward Monte-Carlo equation of motion leads to improved time convergence.

This research was sponsored by the Oak Ridge National Laboratory,managed by UT–Battelle, LLC, for the U.S. Department of Energy undercontract DE–AC05–00OR22725. Research was supported in part by anappointment to the ORNL Postdoctoral Research Associates Program,administered jointly by Oak Ridge National Laboratory and the Oak Ridge Institute for Science and Education.
The author wishes to thank his colleagues in the Fusion Energy Division Radiofrequency Theory Group (Don Batchelor, Lee Berry, Mark Carter, and Fred Jaeger) for helpful comments during the work on this article.

## Appendix A A new perspective on the forward Monte-Carlo method

The aim of this appendix is to derive the conventional, forward, weighted Monte-Carlo method in the same manner as the backward method was derived . The starting point is the *forward* SDE:

$$dX^{}=\mu dt+\sigma dW.$$
(8)
Solving this SDE is trivial; following exactly the same procedure as for the backward SDE, we get

$$X^{}(t+\mathrm{\Delta }t)=X^{}(t)+\mu \mathrm{\Delta }t+\zeta \sigma \sqrt{\mathrm{\Delta }t}+𝒪(\mathrm{\Delta }t),$$
(9)
where $`\zeta `$ is a zero-mean, unit-variance Gaussian random number, $`\zeta N(0,1)`$. When the SDE is backward, the obvious initial condition is $`X^{}(t=T)=x`$. For the forward SDE (8), we simply do not know which initial condition to impose; we will have to leave $`X^{}(t=0)`$ undefined for now.

As before , we use the Itô formula to differentiate $`f(X(t),t)`$ and obtain

$$f(X^{}(T),T)=f(X^{}(0),0)+_0^T2f_t𝑑t+_0^T\sigma f_x𝑑W,$$
(10)
where again we have identified $`\mu =D/x`$, $`\sigma =\sqrt{2D}`$, and used$`f_t+\mu f_x+\frac{1}{2}\sigma ^2f_{xx}=`$ \[Eq. (1)\] $`=2f_t`$. Using the macroscopic initial condition, $`f(x,0)=\mathrm{\Phi }(x)`$, we get

$$f(X^{}(T),T)=\mathrm{\Phi }(X^{}(0))_0^T\sigma f_x𝑑W.$$
(11)
At this stage of the derivation of the *backward* method, we used the microscopic initial condition and found that the equivalent of the LHS of Eq. (11) was in fact a solution to Eq. (1): $`f(X^{}(t=T),T)=`$ \[$`X^{}(t=T)=x`$\] $`=f(x,T)`$. But here, in the forward derivation, we cannot do that because $`X^{}(T)x`$! Going forward in time, $`X^{}(T)`$ is an unknown; we could use the Monte-Carlo equation of motion, Eq. (9), to find $`X^{}(T)`$, but we do not have an initial condition! So, is there any way to get $`f(x,T)`$ from Eq.(8)? The answer is yes; there is a (rather contrived) way.

If we let $`f(x,T)`$ be a distribution, we can write it as

$$f(x,T)=E[f(X^{}(T),T)\times \delta (xX^{}(T))].$$
(12)
Substituting Eq. (11) into Eq. (12), we get:

$$f(x,T)=E[\mathrm{\Phi }(X^{}(0))\times \delta (xX^{}(T))].$$
(13)
The numerical approximation of this expectation value is the forward weighted Monte-Carlo solution:

$$f_{}(x,T)=N^1\underset{i=1}{\overset{N}{}}\mathrm{\Phi }(X_i^{}(0))\times \delta (xX_i^{}(T)).$$
(14)
Now, we still need to know what the microscopic initial condition on $`X_i^{}(0)`$ should be. The best we can do is to make sure that $`f_{}(x,T=0)\mathrm{\Phi }(x)`$. Note that “approximately equal” must be given a very liberal definition because $`f_{}(x,T=0)`$ is a jagged sum of $`\delta `$-functions, whereas $`\mathrm{\Phi }(x)`$ is in general a smooth function. The stochastic variables $`X_i^{}(0)`$ can now be found by following the stochastic trajectories given by the forward Monte-Carlo equation of motion, Eq. (9), with the “fuzzy initial condition”:

$$X_i^{}(0)\{X_i^{}(0)|f_{}(x,T=0)\mathrm{\Phi }(x)\}.$$
(15)
## Appendix B A higher-order Monte-Carlo equation of motion

The diffusive term $`\zeta \sigma (X(t))\sqrt{\mathrm{\Delta }t}`$ of the Monte-Carlo equations of motion, Eqs. (3) and (5), is an approximation of the Itô integral

$$_t^{t+\mathrm{\Delta }t}\sigma (X(s))𝑑W(s)=\sigma (X(t))_t^{t+\mathrm{\Delta }t}𝑑W(s)+𝒪(\mathrm{\Delta }t),$$
(16)
where $`X(s)`$ is the stochastic process that solves the SDE

$$dX=\mu (X(s))ds+\sigma (X(s))dW(s),$$
(17)
where $`W`$ is a Wiener process, and we impose the initial condition $`X(t)=x_t`$. The approximation Eq. (16) replaces $`\sigma (x)`$ with its zero-order Taylor expansion, $`\sigma _0(x)=\sigma (x_t)`$. To reduce the error in Eq. (5) to $`𝒪(\mathrm{\Delta }t^{3/2})`$ we Taylor expand $`\sigma (x)`$ to first order,

$$\sigma _1(x)=\sigma (x_t)+\sigma {}_{}{}^{}(x_t)(xx_t),$$
(18)
and approximate the solution to Eq. (17) with

$$X(s)=x_t+\mu (x_t)(st)+\sigma (x_t)[W(s)W(t)].$$
(19)
The Itô integral over $`\sigma `$ becomes:

$`{\displaystyle _t^{t+\mathrm{\Delta }t}}`$ $`\sigma `$ $`(X(s))dW(s)=\sigma (x_t){\displaystyle _t^{t+\mathrm{\Delta }t}}𝑑W(s)+`$
$`\sigma `$ $`{}_{}{}^{}(x_t){\displaystyle _t^{t+\mathrm{\Delta }t}}\{\mu (x_t)(st)+\sigma (x_t)[W(s)W(t)]\}dW(s)+𝒪(\mathrm{\Delta }t^{3/2})=`$
$`\sigma `$ $`(x_t){\displaystyle _t^{t+\mathrm{\Delta }t}}dW(s)+\sigma (x_t)\sigma {}_{}{}^{}(x_t){\displaystyle _t^{t+\mathrm{\Delta }t}}[W(s)W(t)]dW(s)+𝒪(\mathrm{\Delta }t^{3/2}),`$
where $`\mu (x_t)_t^{t+\mathrm{\Delta }t}(st)𝑑W(s)`$ is $`𝒪(\mathrm{\Delta }t^{3/2})`$ and hence was neglected. The first term on the RHS of Eq. (LABEL:eq:sigma\_integral\_high) is just the low-order approximation Eq. (16). The second term is the next-order correction. To calculate this correction term, we first need some intermediate results.

Following Björk , we introduce the stochastic variables

$$I_n=\underset{j=0}{\overset{n1}{}}W(t_k)[W(t_{k+1})W(t_k)],$$
(21)
and

$$B_n=\underset{j=0}{\overset{n1}{}}W(t_{k+1})[W(t_{k+1})W(t_k)],$$
(22)
where $`t_j=jt/n,j=0,1,\mathrm{},n1`$. It trivially follows that

$$B_n+I_n=\underset{j=0}{\overset{n1}{}}[W(t_{k+1})+W(t_k)][W(t_{k+1})W(t_k)]=W^2(t),$$
(23)
and

$$B_nI_n=\underset{j=0}{\overset{n1}{}}[W(t_{k+1})W(t_k)]^2=S_n(t).$$
(24)
In Appendix A of Ref.  we showed that

$$\underset{n\mathrm{}}{lim}S_n(t)=t.$$
(25)
Subtracting Eq. (24) from Eq. (23) and taking the limit $`n\mathrm{}`$, we obtain

$$\underset{n\mathrm{}}{lim}I_n(t)=_0^tW(s)𝑑W(s)=\frac{1}{2}[W^2(t)t].$$
(26)
Applying Eq. (26) to Eq. (LABEL:eq:sigma\_integral\_high) we get

$`{\displaystyle _t^{t+\mathrm{\Delta }t}}`$ $`\sigma `$ $`(X(s))dW(s)=\sigma (X(t))[W(t+\mathrm{\Delta }t)W(t)]+`$
$`\frac{1}{2}`$ $`\sigma (X(t))\sigma {}_{}{}^{}(X(t))\{[W(t+\mathrm{\Delta }t)W(t)]^2\mathrm{\Delta }t\}+𝒪(\mathrm{\Delta }t^{3/2}).`$
In the absence of advection, as is the case for Eq. (1), $`\sigma \sigma {}_{}{}^{}=\mu `$, and the higher-order Monte-Carlo equations of motion become

$`X_i^{}(t\pm \mathrm{\Delta }t)`$ $`=`$ $`X_i^{}(t)+\mu \mathrm{\Delta }t+\zeta \sigma \sqrt{\mathrm{\Delta }t}+\frac{1}{2}(\zeta ^21)\mu \mathrm{\Delta }t+𝒪(\mathrm{\Delta }t^{3/2})=`$ (28)
$`X_i^{}(t)+\frac{1}{2}(1+\zeta ^2)\mu \mathrm{\Delta }t+\zeta \sigma \sqrt{\mathrm{\Delta }t}+𝒪(\mathrm{\Delta }t^{3/2}),`$
where we have used the fact that $`W(t+\mathrm{\Delta }t)W(t)N(0,\sqrt{\mathrm{\Delta }t})`$; $`\zeta `$ is a zero-mean, unit-variance Gaussian random number; $`\zeta N(0,1)`$; and $`\mu =\mu (X(t))`$ and $`\sigma =\sigma (X(t))`$.

In Ref.  Milshtein derived Eq. (LABEL:eq:sigma\_integral\_high\_eval) by introducing

$$X^{}(t+\mathrm{\Delta }t)=x_t+c_1\mathrm{\Delta }t+c_2[W(t+\mathrm{\Delta }t)W(t)]+c_3[W(t+\mathrm{\Delta }t)W(t)]^2,$$

and finding the $`c_1`$, $`c_2`$ and $`c_3`$ that minimize $`E[(X^{}(t+\mathrm{\Delta }t)X(t+\mathrm{\Delta }t))^2]`$.


format: text

text:
# Statistical Model of Superconductivity in a 2D Binary Boson-Fermion Mixture

## 1 Introduction

Recent experiments  indicate that composite bosons in ultra-cold clouds of most alkali atoms do indeed Bose-Einstein (BE) condense. Since Cooper pairs (CPs) of fermions (electrons or holes) in a many-fermion system form composite bosons in the sense of coupling to integer angular momentum, it is natural to consider the possible BE condensation of such pairs. The belief that some such condensate is central to superconductivity is more than 50 years old . High-$`T_c`$, as well as some organic, superconductors  are quasi-two-dimensional (2D). Quasi-1D superconductors have also been found . BE condensation (BEC) is impossible in two or less space dimensions  for usual or “ordinary” bosons (i.e., with a quadratic energy-momentum, or dispersion, relation). It is however still possible to have BEC in all dimensions $`d>1`$ for non-interacting bosons if they obey a linear dispersion relation —such as CPs moving in the Fermi sea. This possibility arises because the Hohenberg theorem , which prohibits BEC in 2D, relies on an $`f`$-sum rule based on the quadratic dispersion relation appropriate to bosons  moving in a vacuum. Such a linear dispersion relation for the CPs in a binary boson-fermion mixture model was recently found  to be consistent, without any adjustable parameters, with the anomalous linear (quadratic) temperature-dependence above $`T_c`$ in the resistivity of optimally-doped (overdoped) cuprates whether hole- or electron-doped. For the observed quadratic $`T`$-dependence in overdoped samples linear-dispersion CP charge carriers are essential.

Although extensive studies in the BCS-Bose “crossover” problem in superconductivity have spanned  a period of over thirty years, we note that BEC is distinct from the standard (i.e., zero center-of-mass momentum CPs) BCS theory condensation where only that one bosonic state exists.

In this paper it is shown that in addition BEC is still possible in 2D even if the number of composite bosons (pairs of fermions) in a binary boson-fermion mixture is not fixed—as chemical/thermal equilibrium renders it coupling- and temperature-dependent—as long as the total number of fermions is fixed. This gives rise to an interesting statistical-mechanics problem irrespective of the particular mechanism for pair formation, and may have a vital application for superconductivity as well as for (neutral-atom) superfluidity such as in liquid 3He , dilute mixtures of 3He in 4He , or in trapped Fermi gases . The statistical model dealt with here may be seamlessly linked to BCS theory, via the fermionic energy gap, when boson/unpaired-fermion interactions are included as, e.g., in Refs.  and . However, in these two papers the quadratic CP dispersion relation has been assumed. The quadratic form has recently been shown  to apply only in the zero-density or vacuum limit when the Fermi sea disappears.

In Section 2 we recall a 2D gas of fermions at $`T=0`$ interacting via a constant pairing interaction in an annulus about the Fermi surface—viz., the BCS model interaction. The binding energy of a single pair near the Fermi surface (CP problem) decreases practically linearly with the center-of-mass momentum (CMM) of the pair for all values of the momentum below breakup, the breakup momentum typically being only about four orders of magnitude smaller than the Fermi momentum. In Section 3 we discuss why the interacting many-fermion system can be treated as a set of independent CPs (i.e., composite bosons with fermion number two) mixed in with pairable fermions which are not bound into pairs, i.e., unpaired fermions. In Section 4 the more realistic scenario is considered of the BEC of these pairs, incorporating pair breakup beyond a certain CMM. Although the number of pairs is not fixed but rather strongly coupling- and temperature-dependent, BEC is still possible in 2D. A simple binary boson- fermion statistical model is introduced by constructing the Helmholtz free energy for an ideal mixture of pairable but unpaired fermions plus paired fermions (both zero and nonzero CMM pairs), all in chemical and thermal equilibrium. The latter results through extrema of the free energy in: a) the pairable fermion occupation probabilities; b) the excited boson numbers (nonzero CMM CPs); and c) the ground boson number (zero CMM pairs). In Section 5 the coupling- and temperature-dependence of the boson number is derived. In Section 6 the critical BEC singularity temperature is obtained first by ignoring the unpaired fermions in a pure boson-gas model and then exactly for the boson-fermion binary mixture model from a $`T`$-dependent dispersion relation derived and calculated numerically, and results compared with empirical data. Finally, Section 7 gives conclusions.

## 2 Cooper-pair dispersion relation

Consider a 2D system of N fermions of mass $`m`$ confined in a square “pen” of area $`L^2`$ and interacting pairwise via the BCS model interaction

$$V_{𝐤,𝐤^{}}=\{\begin{array}{cc}V& \text{if}\mu (T)\mathrm{}\omega _D<ϵ_{k_1}(\mathrm{}^2k_1^2/2m),ϵ_{k_2}<\mu (T)+\mathrm{}\omega _D\hfill \\ 0& \text{otherwise},\hfill \end{array}$$
(1)
where $`𝐤\frac{1}{2}(𝐤_1𝐤_2)`$ is the relative wavevector of the two particles (see Fig. 3 below); $`V_{𝐤,𝐤^{}}`$ the 2D double Fourier integral of the underlying non-local interaction $`V(𝐫,𝐫^{})`$ in the relative coordinate $`𝐫=𝐫_1𝐫_2;\mu (T)`$ the ideal Fermi gas chemical potential which at $`T=0`$ becomes the Fermi energy $`E_F\mathrm{}^2k_F^2/2m`$ with $`k_F`$ the Fermi wavenumber; $`\mathrm{}\omega _D\mathrm{}^2k_D^2/2m`$ the width of the annulus about the Fermi circle in which the pairing interaction is nonzero, with $`\omega _D`$ being the Debye frequency. This model interaction mimics the net effect of an attractive electron-phonon interaction overwhelming the repulsive interfermion Coulomb repulsions, whenever $`V>0`$.

If $`\mathrm{}𝐊=\mathrm{}(𝐤_1+𝐤_2)`$ is the center-of-mass momentum (CMM) of a pair, let $`E_K`$ be its total energy (besides the CP rest-mass energy). The eigenvalue (CP ) equation for a pair of fermions at $`T=0`$ immersed in a background of $`N2`$ inert, spectator fermions within a (sharp) Fermi circular perimeter of radius $`k_F`$ is then

$$1=V\underset{𝐤}{}^{^{}}\frac{\theta (k_1k_F)\theta (k_2k_F)}{2ϵ_k(E_K\mathrm{}^2K{}_{}{}^{2}/4m)},$$
(2)
where again $`ϵ_k\mathrm{}^2k^2/2m`$, $`\theta (x)`$ is the Heaviside unit step function, and the prime on the summation sign denotes the conditions

$$k_1|𝐤+\frac{1}{2}𝐊|<(k_F^2+k_D^2)^{1/2}\text{and}k_2|𝐤\frac{1}{2}𝐊|<(k_F^2+k_D^2)^{1/2}$$
(3)
ensuring that the pair of fermions above the Fermi “surface” cease interacting beyond the annulus of energy thickness $`\mathrm{}\omega _D`$ in accordance with (1), thereby restricting the summation over $`𝐤`$ for a given fixed $`𝐊`$. Without these restrictions (2) would just be the Schrödinger equation in momentum space for the pair. Setting $`E_K2E_F\mathrm{\Delta }_K`$, the pair is bound if $`\mathrm{\Delta }_K>0`$, and (2) becomes an eigenvalue equation for the (positive) pair binding energy $`\mathrm{\Delta }_K`$. Our $`\mathrm{\Delta }_K`$ and $`\mathrm{\Delta }_0`$ should not be confused with the BCS energy gap $`\mathrm{\Delta }(T)`$.

Let $`\lambda g(E_F)V0`$ be a dimensionless coupling constant with $`g(E_F)`$ the electronic density-of-states (for each spin) at the Fermi surface in the normal (i.e., interactionless) state, which in 2D is constant

$$g(ϵ)=L^2m/2\pi \mathrm{}^2g.$$
(4)
The Cooper equation (2) for the unknown quantity $`\mathrm{\Delta }_K`$ is analyzed in Ref. . For zero CMM, $`K=0`$, it becomes a single elementary integral, with the familiar  solution

$$\mathrm{\Delta }_0=\frac{2\mathrm{}\omega _D}{e^{2/\lambda }1}$$
(5)
valid for all coupling $`\lambda `$. For small $`K`$, it is not too difficult to extract  the asymptotic result

$$\mathrm{\Delta }_K\underset{K0}{}\mathrm{\Delta }_0\frac{2}{\pi }\left[1+\frac{\mathrm{\Delta }_0}{2\mathrm{}\omega _D}(1+\sqrt{1+\nu })\right]\mathrm{}v_FK+O(K^2)\underset{\lambda 0}{}\mathrm{\Delta }_0\frac{2}{\pi }\mathrm{}v_FK+O(K^2)$$
(6)
where $`\nu \mathrm{\Theta }_D/T_F`$, and $`v_F`$ is the Fermi velocity defined through $`E_F\mathrm{}^2k_F^2/2m=\frac{1}{2}mv_F^2`$. For weak coupling, $`\lambda 0`$, this linear dispersion relation gives the 2D analog of the 3D result stated as far back as 1964 in Ref. , p. 33 (see also, Ref. , p. 336) but with the 2D coefficient $`2/\pi `$ of the last expression of (6) replaced by $`1/2`$.

## 3 Justification of boson formalism

These CP boson-like structures could be called “quasi-bosons” since their creation and annihilation operators are known not to obey the usual boson commutation relations , p. 38. However, they do obey the Bose-Einstein distribution since the energy $`E_K`$ of the CP is given only by the total CMM, $`K`$, but is independent of the relative momentum $`k`$. Thus, the possible energy states for the pair are $`E_K`$ as defined in (2). The number of pairs $`N_𝐊`$ that can occupy such a state can take on indefinite values since there exist also indefinitely many relative momenta, namely

$$N_𝐊\underset{𝐤}{}𝒩_{𝐤,𝐊}=0,1,2,\mathrm{}.$$
(7)
Here, $`𝒩_{𝐤,𝐊}=0,1`$ is the number of pairs characterized by both k and K, and is the same number as that characterized by definite $`𝐤_1`$ and $`𝐤_2`$, namely $`𝒩_{𝐤,𝐊}=n_{𝐤_1}n_{𝐤_2}=0,1`$ where $`n_{𝐤_i}=0,1`$ is the occupation number for a single fermion, these remarks all referring to singlet pairing. Much of all this has been known  at least since 1958, albeit in somewhat different language.

This view of an actual Cooper pair should not be confused with, say, an Anderson  phonon-like collective excitation (or modes) with weak-coupling dispersion relation—in 2D  given by $`(1/\sqrt{2})\mathrm{}v_FK`$ in the long-wavelength limit, and which evolves into the plasmon when Coulomb repulsions between fermions are switched on. CPs here, like deuterons, carry fermion number two and as such are definite in number (although in the CP case this number is coupling- and temperature-dependent) and can thus undergo BEC. This is distinct from collective excitations which are indefinite in number. Park , e.g., distinguishes between “permanent” and “ephemeral” bosons, the latter sometimes being referred to as “quasiparticles” to distinguish from the former “particles”.

For $`N_B`$ ordinary bosons of mass $`m_B`$ and energy $`\epsilon _K=\mathrm{}^2K^2/2m_B`$ in any positive dimension, $`d>0`$, a temperature singularity $`T_c`$  appears in the number equation $`N_B=_𝐊[e^{(\epsilon _K\mu _B)/k_BT}1]^1`$ at vanishing bosonic chemical potential $`\mu _B\underset{}{<}0`$ when the number of $`𝐊=0`$ bosons just ceases to be negligible upon cooling. It is given by

$$T_c=\frac{2\pi \mathrm{}^2}{m_Bk_B}\left[\frac{n_B}{g_{d/2}(1)}\right]^{2/d}$$
(8)
with $`n_B`$ the boson particle density $`N_B/L^d`$, and $`g_{d/2}(z)`$ the usual Bose integrals

$$g_\sigma (z)\frac{1}{\mathrm{\Gamma }(\sigma )}_0^{\mathrm{}}𝑑x\frac{x^{\sigma 1}}{z^1e^x1}=\underset{l=1}{\overset{\mathrm{}}{}}\frac{z^l}{l^\sigma }\underset{z1}{}\zeta (\sigma ),$$
(9)
where $`\mathrm{\Gamma }(\sigma )`$ is the gamma function and $`\zeta (\sigma )`$ the Riemann zeta function of order $`\sigma `$. The last identification in (9) holds when $`\sigma >1`$ for which $`\zeta (\sigma )<\mathrm{}`$, while the series $`g_\sigma (1)`$ diverges for $`\sigma 1`$, thus giving $`T_c=0`$ for $`d2`$. For $`d=3`$ one has $`\zeta (3/2)2.612`$ so that (8) becomes the familiar formula $`T_c3.31\mathrm{}^2n_B^{2/3}/m_Bk_B`$ of “ordinary” BEC. On the other hand, for bosons with (positive) excitation energy $`\epsilon _K\mathrm{\Delta }_0\mathrm{\Delta }_K`$ given approximately by the linear term in (6) for all $`K`$, the singularity that lead to (8) now yields , for weak coupling,

$$T_c=\frac{a(d)\mathrm{}v_F}{k_B}\left[\frac{\pi ^{\frac{d+1}{2}}n_B}{\mathrm{\Gamma }(\frac{d+1}{2})g_d(1)}\right]^{1/d}$$
(10)
where  $`a(d)=1,2/\pi `$ and $`1/2`$ for $`d=1,2`$ and $`3`$, respectively. Note that now $`T_c>0`$ for all $`d>1`$, which is precisely the dimensionality range of all known superconductors including the quasi-1D organo-metallic (Bechgaard) salts . This is not inconsistent with the Hohenberg theorem  that there is no broken symmetry, i.e., long-range order, in a Bose fluid for $`d`$ = 1 or 2, since this is based on an $`f`$-sum rule for bosons with a quadratic dispersion relation. Indeed, both (8) and (10) are special cases of of the more general expression  for any space dimensionality $`d>0`$ and any boson dispersion relation $`\epsilon _K=C_sK^s`$ with $`s`$ $`>0`$ and $`C_s`$ a constant, given by

$$T_c=\frac{C_s}{k_B}\left[\frac{s\mathrm{\Gamma }(d/2)(2\pi )^dn_B}{2\pi ^{d/2}\mathrm{\Gamma }(d/s)g_{d/s}(1)}\right]^{s/d}.$$
(11)
In what follows the number of bosons will be temperature-dependent and it is in conserving the fermion number that the singularity arises. As is the case for the pure boson gas, a linear rather than a quadratic dispersion relation will be needed to obtain BEC in 2D. This emerges in a statistical model of an ideal binary mixture of bosons (the CPs) and unpaired (both pairable and unpairable) fermions in chemical equilibrium , for which thermal pair-breaking into unpaired pairable fermions is explicitly allowed.

## 4 First-principles statistical model

Under interaction (1) at any $`T`$ the total number of fermions in 2D is $`N=L^2k_F^2/2\pi =N_1+N_2`$ and is just the number of non-interacting (i.e., unpairable) fermions $`N_1`$ plus the number of pairable ones $`N_2`$. The unpairable fermions obey the usual Fermi-Dirac distribution with fermionic chemical potential $`\mu `$. On the other hand, the $`N_2`$ pairable fermions are simply those in the interaction shell of energy width $`\mathrm{}\omega _D`$ so that

$$N_2=2_{\mu \mathrm{}\omega _D}^{\mu +\mathrm{}\omega _D}𝑑ϵ\frac{g(ϵ)}{e^{\beta (ϵ\mu )}+1}=2g\mathrm{}\omega _D,$$
(12)
since the density of electronic states (4) is constant and the remaining integral exact. At any interfermionic coupling and temperature these fermions form an ideal mixture of pairable but unpaired fermions plus CPs that are created near the single-fermion energy $`\mu (T)`$, with binding energy $`\mathrm{\Delta }_K(T)`$ $`0`$ and total energy

$$E_K(T)2\mu (T)\mathrm{\Delta }_K(T).$$
(13)
This generalizes the $`T=0`$ equation $`E_K2E_F\mathrm{\Delta }_K`$ introduced below (3).

The Helmholtz free energy $`F=ETS`$, where $`E`$ is the internal energy and $`S`$ the entropy, for this binary “composite boson/pairable-but-unpaired-fermion system” at temperatures $`TT_c`$ is then

$`F_2`$ $`=`$ $`2{\displaystyle _{\mu \mathrm{}\omega _D}^{\mu +\mathrm{}\omega _D}}𝑑ϵg(ϵ)\left\{n_2(ϵ)ϵ+k_BT\left[n_2(ϵ)\mathrm{ln}n_2(ϵ)+\{1n_2(ϵ)\}\mathrm{ln}\{1n_2(ϵ)\}\right]\right\}`$ (14)
$`+[2\mu (T)\mathrm{\Delta }_0(T)]N_{B,0}(T)`$
$`+{\displaystyle \underset{K>0}{\overset{K_0}{}}}\{[2\mu (T)\mathrm{\Delta }_K(T)]N_{B,K}(T)`$
$`+k_BT[N_{B,K}(T)\mathrm{ln}N_{B,K}(T)\{1+N_{B,K}(T)\}\mathrm{ln}\{1+N_{B,K}(T)\}]\}.`$
The integral term is the contribution from the unpaired fermions and runs over all levels in the energy shell where the BCS model interaction is nonzero, $`n_2(ϵ)`$ being the average number of unpaired but pairable fermions with energy $`ϵ`$; the prefactor two comes from the spin. The second term gives the free energy of the bosons with CMM $`K=0`$ since their entropy is negligible in the thermodynamic limit; here $`N_{B,0}(T)`$ is the number of (bosonic) CPs with zero CMM at temperature $`T`$. The summation term represents the free energy of the bosons with nonzero CMM, while $`N_{B,K}(T)`$ is that with arbitrary nonzero CMM $`K,`$ and the cutoff $`K_0`$ is defined  by $`\mathrm{\Delta }_{K_0}0`$. The free energy $`F_2`$ is to be minimized subject to the constraint that the total number of pairable fermions $`N_2`$ is conserved.

If $`N_{20}(T)`$ is the number of pairable but unpaired fermions, the relevant number equation for the pairable (i.e., active) fermions is then

$$N_2=N_{20}(T)+2[N_{B,0}(T)+N_{B,0<K<K_0}(T)]N_{20}(T)+2N_B(T),$$
(15)
where $`N_{B,0<K<K_0}(T)`$ denotes the total number of “excited” bosonic pairs (namely with CMM such that $`0<K<K_0`$), i.e., $`N_{B,0<K<K_0}(T)_{0<K<K_0}N_{B,K}(T)`$. Minimizing the free energy, subject to the constraint that (15) be a constant, is equivalent to minimizing the grand potential

$$\mathrm{\Omega }_2=F_2\mu _2N_2.$$
(16)
a) Minimizing $`\mathrm{\Omega }_2`$ with respect to the fermion occupation probabilities $`n_2(ϵ)`$ yields the Fermi-Dirac distribution with fermion chemical potential $`\mu _2`$, not $`\mu `$, namely

$$n_2(ϵ)=\frac{1}{e^{\beta (ϵ\mu _2)}+1};\beta (k_BT)^1.$$
(17)
Thus the total number of pairable (but unpaired) fermions then becomes

$$N_{20}(T)2_{\mu \mathrm{}\omega _D}^{\mu +\mathrm{}\omega _D}𝑑ϵg(ϵ)n_2(ϵ)=2_{\mu \mathrm{}\omega _D}^{\mu +\mathrm{}\omega _D}𝑑ϵ\frac{g(ϵ)}{e^{\beta (ϵ\mu _2)}+1},$$
(18)
and should be compared with (12) for $`N_2`$ which contains only $`\mu `$. Since in 2D $`g(ϵ)`$ is a constant (4), (18) becomes the exact expression

$$N_{20}(T)=\frac{2g}{\beta }\mathrm{ln}\left[\frac{1+e^{\beta (\mu \mu _2\mathrm{}\omega _D)}}{1+e^{\beta (\mu \mu _2+\mathrm{}\omega _D)}}\right].$$
(19)
b) Minimizing $`\mathrm{\Omega }_2`$ with respect to the excited boson numbers $`N_{B,K}(T)`$, $`K>0`$, yields the Bose-Einstein distribution summed over all $`0<K<K_0`$, namely

$$N_{B,0<K<K_0}(T)\underset{K>0}{\overset{K_0}{}}N_{B,K}(T)=\underset{K>0}{\overset{K_0}{}}[e^{\beta \{E_K(T)2\mu _2\}}1]^1.$$
(20)
The factor multiplying $`\beta `$ in (20) may be rewritten as $`\epsilon _K(T)\mu _B(T)`$, where $`\epsilon _K(T)\mathrm{\Delta }_0(T)\mathrm{\Delta }_K(T)0`$ is a (nonnegative) excitation energy as suggested by (6), while $`\mu _B(T)`$ turns out to be

$$\mu _B(T)=2\left[\mu _2(T)\mu (T)\right]+\mathrm{\Delta }_0(T).$$
(21)
This allows rewriting (20) in the more meaningful boson form

$$N_{B,0<K<K_0}(T)=\underset{K>0}{\overset{K_0}{}}[e^{\beta \{\epsilon _K(T)\mu _B(T)\}}1]^1$$
(22)
where $`\mu _B(T)`$ is clearly the bosonic chemical potential associated with the entire binary mixture.

c) Finally, minimizing $`\mathrm{\Omega }_2`$ with respect to the number of zero CMM (or, “ground state”) bosons $`N_{B,0}(T)`$ gives

$$2[\mu _2(T)\mu (T)]+\mathrm{\Delta }_0(T)=0(0TT_c),$$
(23)
valid only in the stated temperature range as $`N_{B,0}(T)`$ is negligible for all $`T>T_c`$. However, in view of (21) this implies that $`\mu _B(T)`$ $`=0`$ for all $`0TT_c`$—which is precisely the BEC condition for a pure boson gas, even though one now deals with a binary boson-fermion mixture.

## 5 Boson number

To determine $`N_B(T)`$ from (15) we need (19) which with (23) reduces to

$$N_{20}(T)=\frac{2g}{\beta }\mathrm{ln}\left[\frac{1+e^{\beta \{\mathrm{\Delta }_0(T)/2\mathrm{}\omega _D\}}}{1+e^{\beta \{\mathrm{\Delta }_0(T)/2+\mathrm{}\omega _D\}}}\right](0TT_c).$$
(24)
At $`T=0`$ two distinct coupling regimes emerge by inspecting (24): a) for $`\mathrm{\Delta }_0/2<\mathrm{}\omega _D`$ or, from (5) for $`\lambda 2/\mathrm{ln}22.89`$, we have that $`N_{20}(0)=2g(\mu )(\mathrm{}\omega _D\mathrm{\Delta }_0/2)`$; while b) for $`\mathrm{\Delta }_0/2>\mathrm{}\omega _D`$ (or $`\lambda 2.89`$) $`N_{20}(0)`$ is identically zero. Hence, the number of bosons $`N_B(0)`$ at $`T=0`$ from (15) is just $`N_B(0)=\frac{1}{2}[N_2N_{20}(0)]`$. Using (12) for $`N_2`$ the fractional number of pairable fermions that are actually paired at $`T=0`$, namely $`2N_B(0)/N_2=1N_{20}(0)/N_2`$, becomes simply

$$2N_B(0)/N_2=\{\begin{array}{cc}\mathrm{\Delta }_0/2\mathrm{}\omega _D=(e^{2/\lambda }1)^1\underset{\lambda 0}{}e^{2/\lambda }\hfill & (\mathrm{for}\lambda 2/\mathrm{ln}22.89)\hfill \\ 1\hfill & (\mathrm{for}\lambda 2/\mathrm{ln}22.89).\hfill \end{array}$$

This fraction is plotted against coupling $`\lambda `$ in Fig. 1 as $`2n_B(0)/n_2`$, since $`n_B(T)N_B(T)/L^2`$ and $`n_2N_2/L^2`$. Since $`N_B(0)=\frac{1}{2}g\mathrm{\Delta }_0`$ for $`\lambda 2.89`$, only those fermions in an energy shell of width $`\mathrm{\Delta }_0/2`$ around the Fermi surface actually pair at $`T=0`$, while for $`\lambda 2.89`$ all pairable fermions actually pair up since then $`N_B(0)=g\mathrm{}\omega _D\frac{1}{2}N_2`$. This result contrasts sharply with the “heuristic model” , Eq. (22), where $`2N_B(0)/N_2`$ $`1`$ for all coupling, and is more in line with BCS theory which implies, in any $`d`$, a coupling-dependent fraction estimated (Ref.  p. 128; see also ) to be $`[g(E_F)2\mathrm{\Delta }/2g(E_F)\mathrm{}\omega _D]^2=(\mathrm{\Delta }/\mathrm{}\omega _D)^2(\mathrm{sinh}1/\lambda )^2\underset{\lambda 0}{}4e^{2/\lambda }`$, where $`\mathrm{\Delta }\mathrm{}\omega _D/\mathrm{sinh}(1/\lambda )`$ (again, not to be confused with the CP binding energy $`\mathrm{\Delta }_0`$) is the $`T=0`$ BCS energy gap for the same BCS model interaction (1) used in this paper; this is graphed as the thin curve in Fig. 1 and is seen to be much larger than (5) for fixed $`\lambda `$. The breakdown of BCS theory for BCS model interaction couplings larger than $`\lambda 1.13`$ is clear both because: a) the alluded fraction cannot exceed unity; and b) physically, if the fermionic energy gap $`\mathrm{\Delta }\mathrm{}\omega _D`$ no pairable fermions are available at all. This breakdown is indicated by the dashed curve in Fig. 1. (A strong-coupling many-body model differing from that of BCS theory but based on the BCS model interaction has been solved by Thouless ).

Also displayed in Fig. 1 are two finite-temperature results for $`2N_B(T)/N_2=1N_{20}(T)/N_2`$ which are obtainable from (24) for any $`T`$ provided one knows $`\mathrm{\Delta }_0(T)`$ for any $`T>0`$. For $`T>0`$, the $`\theta (k_1k_F)\theta (ϵ_{k_1}E_F)`$ in (2) becomes $`1n(\xi _{k_1})`$, where $`n(\xi _{k_1})(e^{\beta \xi _{k_1}}+1)^1`$ is the Fermi-Dirac distribution with $`\xi _{k_1}ϵ_{k_1}\mu (T)`$, with the ideal fermion gas chemical potential $`\mu (T)`$ in 2D being given exactly by

$$\mu (T)=\beta ^1\mathrm{ln}(e^{\beta E_F}1)\underset{T0}{}E_F.$$
(25)
Note that $`\mu (T)`$ decreases monotonically with temperature from its maximum value of $`E_F`$ but does not turn negative until $`T=T_F/\mathrm{ln}21.44T_F`$ so that the BCS model interaction (1), which requires $`\mu (T)`$ to be nonnegative, will not break down (i.e., become meaningless) over the entire range of temperatures relevant in this paper, see Fig. 4 below. Similar arguments hold for $`\theta (k_2k_F)`$. Since $`k_1=k_2`$ implies that $`\xi _{k_1}=\xi _{k_2}`$, (2) then leads to a simple generalization to finite-temperature of the $`K=0`$ CP equation, namely

$$1=\lambda _0^{\mathrm{}\omega _D}𝑑\xi (e^{\beta \xi }+1)^2[2\xi +\mathrm{\Delta }_0(T)]^1.$$
(26)
Its numerical solution for $`\mathrm{\Delta }_0(T)`$ is illustrated in Fig. 2 for $`\lambda gV=1/2`$ and $`\nu \mathrm{}\omega _D/E_F=0.05`$. Note that if one assumes a $`T^{}`$ such that $`\mathrm{\Delta }_0(T^{})=0`$, the resulting integral in (26) diverges and the equation can only be satisfied for $`\lambda =0`$; thus, there is no temperature $`T^{}`$ at which “depairing” will occur for any fixed $`\nu `$ and any nonzero $`\lambda `$.

## 6 Critical temperature

Neglecting the background unpaired fermions and modeling our system as a pure boson gas of CPs but with temperature-dependent number density $`n_B(T)`$, one converts the explicit $`T_c`$-formula (10) into an implicit one by allowing $`n_B`$ to be $`T`$-dependent. For $`d=2`$ (10) becomes, since $`g_2(1)\zeta (2)=\pi ^2/6`$,

$$T_c=\frac{4\sqrt{3}}{\pi ^{3/2}}\frac{\mathrm{}v_F}{k_B}\sqrt{n_B(T_c)}.$$
(27)
This requires $`n_B(T)N_B(T)/L^2`$ which in turn requires (24), along with $`\mathrm{\Delta }_0(T)`$ as determined from (26), and is given by the expression $`2N_B(T)/N_2=1N_{20}(T)/N_2`$. Solving (27) self-consistently with $`\lambda =1/2`$ gives the remarkably constant value $`T_c/T_F0.004`$ over the entire range of $`\nu \mathrm{}\omega _D/E_F`$ values $`0.030.07`$ typical  of cuprate superconductors. On the other hand, the BCS formula $`T_c^{BCS}1.13\mathrm{\Theta }_De^{1/\lambda }`$ with $`\lambda =1/2`$ gives $`T_c/T_F`$ = 0.005, 0.008 and 0.011 for $`\nu `$ = 0.03, 0.05 and 0.07, respectively. Clearly, both sets of predictions are somewhat small compared with empirical cuprate values of $`T_c/T_F`$ that range  from $`0.010.1`$.

To obtain the exact critical temperature without neglecting the background unpaired fermions, one needs the exact CP excitation energy dispersion relation $`\epsilon _K(T)\mathrm{\Delta }_0(T)\mathrm{\Delta }_K(T)`$ which is neither exactly linear in $`K`$ nor independent of $`T`$. To determine $`\mathrm{\Delta }_K(T)`$ we need a working equation that generalizes Ref.  for $`T>0`$ via the new CP eigenvalue equation (26). Because of symmetry, see Fig. 3, one can restrict the angle $`\theta `$ to the interval $`(0,\pi /2)`$ where $`k_1k_2`$, i.e., to quadrant I. Recalling (13), in $`d`$-dimensions (2) becomes

$$1=V\left(\frac{L}{2\pi }\right)^d^{}𝑑𝐤\frac{\left[1n(\xi _{k_1})\right]\left[1n(\xi _{k_2})\right]}{\mathrm{}^2(k^2k_\mu ^2)/m+\mathrm{\Delta }_K(T)+\mathrm{}^2K^2/4m}.$$
(28)
Here $`k_\mu `$ is such that $`\mu \mathrm{}^2k_\mu ^2/2m`$ and becomes $`k_F`$ as $`T0`$, while $`k_D`$ is such that $`\mathrm{}\omega _D\mathrm{}^2k_D^2/2m`$. The prime on the integral sign now denotes the restrictions

$`k_2^2|𝐤\frac{1}{2}𝐊|^2`$ $`=`$ $`k^2kK\mathrm{cos}\theta +\frac{1}{4}K^2>k_\mu ^2,`$ (29)
$`k_1^2|𝐤+\frac{1}{2}𝐊|^2`$ $`=`$ $`k^2+kK\mathrm{cos}\theta +\frac{1}{4}K^2<k_\mu ^2+k_D^2.`$ (30)
In Fig. 3 the darkest shading corresponds to these (BCS model interaction) restrictions. The conditions (29) and (30) can be studied separately but must be satisfied simultaneously. If $`K<2\sqrt{k_\mu ^2k_D^2}`$, (29) and (30) are equivalent to

$$\begin{array}{ccccc}(k_\mu ^2\frac{1}{4}K^2\mathrm{sin}^2\theta )^{1/2}+\frac{1}{2}K\mathrm{cos}\theta & <& k& <& [(k_\mu ^2+k_D^2)\frac{1}{4}K^2\mathrm{sin}^2\theta ]^{1/2}\frac{1}{2}K\mathrm{cos}\theta .\end{array}$$
(31)
Note that for $`K>\sqrt{k_\mu ^2+k_D^2}\sqrt{k_\mu ^2k_D^2}`$ there exists a minimum value $`\theta _{\mathrm{min}}`$ of $`\theta `$ given by

$$\mathrm{cos}\theta _{\mathrm{min}}\frac{k_D^2}{K\sqrt{2(2k_\mu ^2+k_D^2)K^2}},$$
(32)
while $`\theta _{\mathrm{min}}`$ = 0 for $`K<\sqrt{k_\mu ^2+k_D^2}\sqrt{k_\mu ^2k_D^2}`$. We introduce the dimensionless variables

$$\kappa \frac{K}{2(k_F^2+k_D^2)^{1/2}}1,\xi \frac{k}{k_F},\stackrel{~}{\mathrm{\Delta }}_\kappa \frac{\mathrm{\Delta }_K}{E_F},\nu \frac{\mathrm{\Theta }_D}{T_F}\frac{k_D^2}{k_F^2},$$
(33)
with $`k_B\mathrm{\Theta }_D\mathrm{}\omega _D\mathrm{}^2k_D^2/2m`$ and$`k_BT_FE_F,`$ where $`k_B`$ is Boltzmann’s constant. Recall the $`d=2`$ constant expression (4) for $`g(ϵ)`$, the restrictions (31), and that for $`K0`$ and $`T>0`$ the step functions in (2) $`\theta (k_{1,2}k_F)\theta (|\frac{1}{2}𝐊\pm 𝐤|k_F)`$ become $`[\mathrm{exp}\{\beta [\mathrm{}^2(\frac{1}{2}𝐊\pm 𝐤)^2/2m\mu (T)]\}+1]^1`$—but with $`2ϵ_k`$ in (2) replaced by $`ϵ_{k_1}+ϵ_{k_2}`$, $`E_F`$ by $`\mu (T)`$ and $`\mathrm{\Delta }_K`$ by $`\mathrm{\Delta }_K(T)`$. One finally arrives at a working equation for the binding energy $`\mathrm{\Delta }_K(T)`$ that generalizes Eq. (18) of Ref. , namely

$`1`$ $`=`$ $`{\displaystyle \frac{4}{\pi }}\lambda {\displaystyle _{\theta _{min}}^{\pi /2}}𝑑\theta {\displaystyle _{\xi _{min}(\theta )}^{\xi _{max}(\theta )}}𝑑\xi \xi {\displaystyle \frac{[1+\mathrm{exp}\{\stackrel{~}{\beta }[\xi ^2+(1+\nu )\kappa ^2+2\sqrt{1+\nu }\kappa \xi \mathrm{cos}\theta 1]\}]^1}{2\xi ^2+2(1+\nu )\kappa ^22+\stackrel{~}{\mathrm{\Delta }}_\kappa (\stackrel{~}{T})}}`$ (34)
$`\times [1+\mathrm{exp}\{\stackrel{~}{\beta }[\xi ^2+(1+\nu )\kappa ^22\sqrt{1+\nu }\kappa \xi \mathrm{cos}\theta 1]\}]^1,`$
where $`\nu \mathrm{}\omega _D/\mu `$, $`\xi _{min}(\theta )\sqrt{1+\nu }\kappa \mathrm{cos}\theta +\sqrt{1(1+\nu )\kappa ^2\mathrm{sin}^2\theta }`$, $`\xi _{max}(\theta )`$
$`\sqrt{1+\nu }\kappa \mathrm{cos}\theta +\sqrt{(1+\nu )(1\kappa ^2\mathrm{sin}^2\theta )}`$ and

$$\theta _{min}=\{\begin{array}{cc}& 0\text{if}2\kappa <1\sqrt{(1\nu )/(1+\nu )},\hfill \\ & \mathrm{cos}^1(\nu /\{4\sqrt{1+\nu }\kappa \sqrt{1+\nu /2(1+\nu )\kappa ^2}\})\mathrm{otherwise}.\hfill \end{array}$$

In (34) we have introduced the more general dimensionless quantities $`\xi k/k_\mu `$, $`\stackrel{~}{\mathrm{\Delta }}_\kappa (\stackrel{~}{T})\mathrm{\Delta }_K(T)/\mu `$, where $`\stackrel{~}{T}k_BT/\mu `$ or $`\stackrel{~}{\beta }\mu \beta `$, and $`\kappa K/2\sqrt{k_\mu ^2+k_D^2}`$.

To obtain the critical temperature from the finite-temperature dispersion relation, besides solving (28) for $`\mathrm{\Delta }_K(T)`$, one needs (12), (15), (22) and (25). At $`T=T_c`$ both $`N_{B,0}(T_c)0`$ and $`\mu _B(T_c)0`$ so that one gets the implicit $`T_c`$-equation for the binary mixture gas

$$1=\frac{\stackrel{~}{T}_c}{\nu }\mathrm{ln}\left[\frac{1+e^{\{\stackrel{~}{\mathrm{\Delta }}_0(\stackrel{~}{T}_c)/2\nu \}/\stackrel{~}{T}_c}}{1+e^{\{\stackrel{~}{\mathrm{\Delta }}_0(\stackrel{~}{T}_c)/2+\nu \}/\stackrel{~}{T}_c}}\right]+\frac{8(1+\nu )}{\nu }_0^{\kappa _0(\stackrel{~}{T}_c)}𝑑\kappa \frac{\kappa }{e^{[\stackrel{~}{\mathrm{\Delta }}_0(\stackrel{~}{T}_c)\stackrel{~}{\mathrm{\Delta }}_\kappa (\stackrel{~}{T}_c)]/\stackrel{~}{T}_c}1}.$$
(35)
This must be solved numerically for the exact $`T_c`$ for each $`\lambda `$ and $`\nu `$ in conjunction with (26) for $`\stackrel{~}{\mathrm{\Delta }}_0(\stackrel{~}{T})`$ and (34) for both $`\stackrel{~}{\mathrm{\Delta }}_\kappa (\stackrel{~}{T})`$ and $`\kappa _0(\stackrel{~}{T}_c)`$. Results for $`\lambda =1/2`$ are shown in Table 1 and Fig. 4 for a range of $`\nu `$ values typical  of cuprates.

In order to make this comparison we have taken $`T_\mu /T_F1`$, a very good approximation up to the highest temperatures dealt with. For example, from Fig. 4 the highest $`T_c/T_F0.14`$ already gives $`T_\mu /T_F0.9999`$ from (25), while for smaller $`T_c/T_F`$ the values of $`T_\mu /T_F`$ are even closer to 1. The $`T_c`$ resulting from the exact dispersion relation for $`T=0`$ (dot-dashed curve) is somewhat higher than the exact result (full curve) but lower than that using the linear approximation for $`\mathrm{\Delta }_K(T)`$ (dotted curve). It is also clear that the effect of using the exact or linear (in $`K`$) cases dominates the effect of the dispersion relation $`T`$-dependence. For cuprates $`d2.03`$ has been suggested  to be more realistic as it reflects inter-CuO-layer couplings but our results in that case would be very similar to those reported here for $`d=2`$.

Thus, for $`\nu =0.05`$ the exact $`T_c`$ is seen to be about 46% lower than the heuristic result found in Ref. , Eqs. (15) and (23). It is curious that all results depend very weakly on the $`T`$-dependence of the CP binding energy $`\mathrm{\Delta }_K(T)`$, in spite of its being substantial throughout the temperatures spanned in this paper, as seen in Fig. 2.

We defer study of the condensate fraction $`N_{B,0}(T)/N_B(T)`$ below $`T_c`$ and merely surmise that it may ultimately help explain the apparent absence  in cuprates of the Hebel-Slichter peak of nuclear-spin (NMR) relaxation rates vs temperature for $`0TT_c`$. Such a peak, originally seen  in aluminum, is perhaps the most stringent and qualitatively convincing experimental test of BCS theory (Refs. , p. 71 and , p. 79 ff). Besides cuprates, it is also absent  in several quasi-1D Bechgaard  and in several quasi-2D (ET) organic salt superconductors.

## 7 Conclusions

A simple statistical model treating CPs as non-interacting bosons in thermal and chemical equilibrium with unpaired fermions is proposed. The model gives rise to a boson number that is strongly coupling- and temperature-dependent. Since the CP dispersion relation is approximately linear, it exhibits a Bose-Einstein condensation of zero-CMM pairs at precisely two dimensions. Exact transition temperatures for the boson-fermion mixture based upon the exact CP dispersion relation are in reasonable agreement with empirical cuprate data.

Needless to say, further corrections are yet to be included in the present simple binary mixture boson-fermion model, e.g., i) realistic Fermi surfaces, ii) Van Hove singularities  or other means of accounting for periodic-crystalline effects, as well as iii) the all important $`d`$-wave interfermionic interaction, iv) the boson-fermion interaction and v) residual interbosonic interactions. As to the latter, also generally neglected in BCS theory, if the lowering  of $`T_c`$ in liquid 4He by about 29% with respect to the ideal Bose gas BEC $`T_c`$ is any guide, interbosonic interactions will also lower $`T_c`$ in a more realistic picture. As to the boson-fermion interaction, it is precisely this ingredient that enabled T.D. Lee and coworkers , and Tolmachev  more generally, to link BCS and BEC through a relation stating that the BE condensate fraction is proportional to the (BCS-like) fermionic gap $`\mathrm{\Delta }(T)`$ squared.

## Acknowledgments

We thank A. Salazar for help with Figure 1. M.C., A.P. and A.R. are grateful for partial support from grant PB98-0124, and M.deLl. from grants PB92-1083 and SAB95-0312, both by DGICYT (Spain), and PAPIIT (Mexico) IN102198-9 as well as from CONACyT (Mexico) 27828-E. R.M.Q. and N.J.D. acknowledge the support of the FRD (South Africa). M.deLl. thanks D.M. Eagles, M. Fortes, O. Rojo, and A.A. Valladares for numerous discussions; V.V. Tolmachev for extensive correspondence; A. Polls and R. Escudero for calling his attention to Refs.  and , respectively; and A.N. Kocharian for reading and commenting the manuscript.

Figure Captions

Figure 1. Fractional number of pairable fermions that are actually paired, at three different temperatures, vs. coupling $`\lambda `$ for the present first-principles model (5) (thick curves) and estimated for BCS theory at $`T=0`$ as explained below (5) (thin curve). The number of pairable fermions with the BCS model interaction used is just (12); all of them are actually paired at $`T=0`$ in the heuristic BEC model, Ref.  Eq. (23).

Figure 2. Temperature dependence of $`K=0`$ CP binding energy $`\mathrm{\Delta }_0(T)`$ obtained numerically from (26) for $`\lambda =1/2`$ and $`\nu =0.05`$. Note that when $`T=\mathrm{}`$ (26) is analytical for $`\mathrm{\Delta }_0(\mathrm{})`$; the latter then turns out to be about $`10^8`$, so that the curve saturates from above to this value at $`T=\mathrm{}`$.

Figure 3. Cross-section of overlap “volume” in momentum space (darkest shading) where the tip of the relative wavevector $`𝐤`$ (for two fermions with wavevectors $`𝐤_\mathrm{𝟏}`$ and $`𝐤_\mathrm{𝟐}`$) must point for the attractive BCS model interaction (1) between them to be nonzero and form a Cooper pair of CMM magnitude $`\mathrm{}K`$.

Figure 4. Critical BEC temperature $`T_c`$ in units of $`T_F`$, resulting for the boson-fermion mixture from (35) for $`\lambda =1/2`$ for varying $`\nu \mathrm{}\omega _D/\mu \mathrm{\Theta }_D/T_F`$: with no approximations (full curve); using $`\mathrm{\Delta }_K(T)`$ evaluated at $`T=0`$ (dot-dashed); using the linear-in-$`K`$ approximation for $`\mathrm{\Delta }_K(T)`$ (dotted). The dashed straight line is the BCS formula $`T_c1.13\mathrm{\Theta }_De^{1/2\lambda }`$ for $`\lambda =1/2`$. The very lowest full horizontal line is the solution of the implicit $`T_c`$-equation (27) for the pure unbreakable-boson gas for $`\nu =0.03,0.05`$ and $`0.07`$. A huge $`T_c`$ enhancement thus follows from the mere presence of background unpaired fermions. Cuprate data are taken from Ref. .