Continuous unitary transformation approach to pairing interactions in statistical physics

We apply the flow equation method to the study of the fermion systems with pairing interactions which lead to the BCS instability signalled by the appearance of the off-diagonal order parameter. For this purpose we rederive the continuous Bogoliubov transformation in a fashion of renormalization group procedure where the low and high energy sectors are treated subsequently. We further generalize this procedure to the case of fermions interacting with the discrete boson mode. Andreev-type interactions are responsible for developing a gap in the excitation spectrum. However, the long-range coherence is destroyed due to strong quantum fluctuations.


Introduction
The flow equation method introduced by Wegner [1] and independently proposed by Wilson and G lazek [2] turned out to be a very useful theoretical tool in the condensed matter and high energy physics (for a recent survey of applications see for instance [3]).The main idea behind this technique is to transform the Hamiltonians to diagonal or block-diagonal structure through the continuous unitary transformation upon flexibly adjusting the generating operator.
Such continuous diagonalization can be thought of as a process of renormalizing the coupling constants in a way similar to the Renormalization Group (RG) approach [4].However, the difference, in comparison to the RG scheme, is that instead of integrating out the high energy excitations (fast modes) one renormalizes them in the beginning of the transformation and afterwards the low energy excitations (slow modes) are transformed.We should emphasize that in the quantum mechanics these low energy excitations determine most of the measurable physical properties, the proper description of which is extremely important [5].
In the first parts of this paper we show the way of constructing the continuous unitary transformation for a reduced BCS model.Its exact solution can be obtained using various methods, including the standard single step Bogoliubov transformation [6].Here we rederive this rigorous solution using the continuous unitary method and illustrate how the pairing instability gradually builds up as an asymptotic fixed point.Recently there appeared in the literature several concepts [7,8] where the authors addressed the Cooper pairing by means of the functional RG.Since the continuous unitary transformation belongs to the same family of RG treatments, it would be instructive to adopt Wegner's procedure [9] in order to deal with the symmetry broken states induced by the pairing interactions.
In the second part we generalize this transformation to the case where fermions interact with a single boson mode.Such a situation can be encountered for example when the conduction band electrons are scattered by the Andreev type interaction on the localized electron pairs [10].On some more general grounds it can also refer to the boson type modes coupled with the correlated electrons [11].Another realization of this scenario is possible in ultracold gases where fermion atoms interact with the weakly bound molecules through the Feshbach resonance which finally leads to the atomic superfluidity [12].c T.Doma ński T.Doma ński

The BCS pairing
To have a simple example illustrating how the fast and slow energy modes are scaled during the continuous unitary transformation we first consider the exactly solvable bilinear Hamiltonian which describes the fermions coupled to the pairing field ∆ k .In (1) we use the standard notation for the creation (annihilation) operator ĉ † kσ (ĉ kσ ) and energy ξ k = ε k −µ is measured from the chemical potential µ.Pairing field can be thought of, for instance, as a mean field approximation to the weak pairing potential between electrons ∆ k = − q V k,q ĉ−q↓ ĉq↑ , where V k,q < 0. In general, one can treat ∆ k as the boson operator, e.g., arising from the Hubbard Stratonovich transformation for the interacting fermion system.It can also correspond to a boson whose exchange mediates the pairing between electrons.Coupling to the boson operator will be discussed in the last section.
Hamiltonians of the fermion and/or boson systems having a bilinear structure (1) have been already considered using the flow equation method [1,9,13,14].However, the authors have focused so far on the determination of quasiparticle energies.Here we supplement their analysis by calculating the correlation functions that contain information about the order parameter as well as characteristics of the excitation spectrum.
It is well known that the bilinear Hamiltonian (1) can be diagonalized by the (single step) Bogoliubov transformation [6].All physical quantities, static and dynamic, can be exactly determined.In this work we get the same rigorous solution in a process of continuous diagonalization of the coupled states |k, ↑ and | − k, ↓ .How fast this can be achieved depends on the distance from the Fermi surface |k − k F | (or on the energy ξ k ).Simultaneously with diagonalization there emerge the coherence factors u k , v k .Their final values depend on the relative distance of the momentum k from the Fermi surface.
We construct the continuous unitary transformation Û(l), where l stands for a formal flow parameter varying between the initial l = 0 value and the value l > 0 such that the transformed Hamiltonian Ĥ(l) = Û(l) Ĥ Û−1 (l) is simplified to the needed structure (diagonal, block-diagonal, tridiagonal or any other).With an increase of the flow parameter l, the Hamiltonian evolves according to the general flow equation [1] d Ĥ(l) dl = [η(l), Ĥ(l)], where η(l) ≡ d Û (l) dl Û−1 (l).We now require the transformed Hamiltonian Ĥ(l) to preserve its bilinear structure (1) and let only the coupling constants become l-dependent (renormalized) with the initial conditions ξ k (0) = ξ k , ∆ k (0) = ∆ k and const(0) = 0.This strategy reminds us the scheme of RG approach and it has been shown [1] that the flow parameter l roughly corresponds to 1/ √ Λ [3].There are several ways of choosing the generating operator η(l).In order to reduce the Hamiltonian to a diagonal structure, we follow the Wegner's proposal [1] η(l) = [ Ĥ0 (l), Ĥ(l)], where Ĥ0 (l) = k,σ ξ k (l)ĉ † kσ ĉkσ denotes a diagonal part.With this choice, the Hamiltonian becomes diagonal in the limit l → ∞.An explicit form of the generating operator for (3) is given by η Substituted into the equation (2) it gives which is identical to the following set of flow equations and dconst(l Formally, the solution of the equation ( 7) can be given as which yields that lim l→∞ ∆ k (l) = 0 for all k = k F .When we multiply the equation ( 6) by ξ k (l) and the equation ( 7) by ∆ k (l) * , their sum gives the following invariance This constraint, together with ∆ k (l) vanishing in the limit l → ∞, implies the known Bogoliubov spectrum In figure 1 we plot the pairing field ∆ k (l) and the change (renormalization) of fermion energies ξ k (l) − ξ k at several stages of the continuous transformation.We notice that fast modes (the states distant from the Fermi energy) are transformed in the first part of the process and the change of their energies is rather small.The slow modes (i.e. the low energy excitations) have to be worked on much longer.Asymptotically, at l → ∞, the entire spectrum reduces to the Bogoliubov structure (10).Now we turn our attention to the dynamical quantities which can be expressed via the normal ĉkσ ; ĉ † kσ ω and the anomalous single particle Green's function ĉk↑ ; ĉ−k↓ ω .We introduced the standard Fourier transforms for the retarded fermion Green's function dωe iωt Â; B ω ≡ iΘ(t) Â(t) B + B Â(t) and time evolution is given by Â(t) = e it Ĥ Âe −it Ĥ .The statistical averaging for an arbitrary observable Ô is defined as where β −1 = k B T .Since the trace is invariant under the unitary transformations we can write down where Ô(l) = Û(l) Ô Û−1 (l).It is convenient to compute the trace in the limit l = ∞ when Hamiltonian becomes diagonal.However, the price which we pay for having the simplified Hamiltonian at l = ∞ is the extra necessity to transform the observables Ô → Ô(l) → Ô(∞).This is to be done similarly to the Hamiltonian (2) using With the choice of the generating operator η(l) given by ( 4) we find from the commutator [η(l), ĉk↑ ] that for l > 0 the annihilation operator ĉk↑ gets convoluted with ĉ † −k↓ .Thus, it is natural to propose the following Ansatz for the l-dependent operator Similar analysis for ĉ † −k↓ operator leads to with the initial conditions u k (0) = 1, v k (0) = 0.These parametrized equations (13,14) substituted into (12) with η(l) given in (4) lead to the coupled flow equations After straightforward algebra for (15,16) we obtain the following invariance |u k (l)| 2 + |v k (l)| 2 = const = 1 which assures that the l-dependent fermion operators (13,14) preserve the anticommutation relations Without a loss of generality we assume ∆ k (l) to be real so that the unknown coefficients u k (l), v k (l) become real too.To determine their asymptotic l = ∞ values we can rewrite (15) as and we integrate both sides in the limits 2 we get for the l.h.s.
due to u k (0) = 1.Using equation (7) we can replace 2ξ k (l)∆ k (l)dl by −d∆ k (l)/2ξ k (l) and from the invariance (9) we obtain that The r.h.s. of equation ( 18) gives after integration because ∆ k (∞) = 0. Combining the results (19,20) we obtain and finally we determine the l = ∞ factors In the same way one can show that 2u  In figure 2 we present the coherence factor u 2 k (l) as a function of the flow parameter l.Again we notice that for the momenta distant from the Fermi surface the coherence factors evolve rather quickly from the initial value u 2 k (0) = 1 to the asymptotic values (see the inset).The coherence factors of the slow modes establish later on, similarly to renormalization of the corresponding excitation energies shown in figure 1. Saturation occurs around l ∝ 1/ξ k (∞) 2 .
Since the transformed Hamiltonian Ĥ(∞) is diagonal, we can easily derive the single particle Green's functions (and the higher order Green's functions as well).With the parameterizations (13,14) we obtain ĉk↑ ; ĉ † The equal time expectation values defined by is the Fermi distribution function) yield the following equations for the average occupancy and for the off-diagonal order parameter These equations (25,26) exactly reproduce the rigorous solution of the reduced BCS Hamiltonian.
In the next section we generalize this treatment to account for the scattering of finite momentum fermion pairs.

Finite momentum fermion pairs
The continuous Bogoliubov transformation discussed in the preceding section can be extended to the more general bilinear Hamiltonian (27) eventually reduces to (1) when G k,q = δ q,0 ∆ k .For q = 0, this model allows for scattering of the finite momentum fermion pairs.In what follows below we analyze the effect of such finite momentum scattering.
The generator from Wegner's proposal applied to the flow equation ( 2) induces the off-diagonal terms c † kσ c p =kσ .In order to eliminate them from the transformed Hamiltonian H(l) we replace (28) by After simple algebraic calculations we find that the off-diagonal terms are exactly cancelled if The modified generating operator (29) has a virtue to preserve the initial structure of the Hamiltonian Ĥ We were not able to analytically solve the equations (32)-(34).Therefore, we explored them numerically assuming the tight binding dispersion ξ k = −2t cos(k x a) and using the pairing potential G k,q in a Lorentian form The normalization factor N (n) was taken such that q G k,q = ∆ and ∆/4t = 0.01.
The flow of the operators c which generalizes the previous equations (13,14) due to finite momentum scattering.Substituting (36,37) into the flow equation ( 12) we obtain From these equations ( 38)-(39) we derive the invariance which again leads to the anticommutation relations (17).
Using the parameterization (36,37) we can determine the single particle Green's functions.Let us first consider the diagonal part (in the Nambu representation).The spectral function The usual structure with two delta peaks describes the long-lived quasiparticle modes of energies ±ξ k (∞) similar to the Bogoliubov modes (10) discussed previously for G k,q = ∆ k δ q,0 .The remaining part corresponds to the background spectrum which is continuously spread over a certain energy regime.Due to the invariance (40), the total spectral weight properly obeys the sum rule ∞ −∞ dωA(k, ω) = 1.In figure 3 we plot the local density of states ρ(ω) = k A(k, ω).We notice that the excitation spectrum is no longer truly gaped.There appears a partial suppression of the fermion states around the Fermi energy (corresponding to ω = 0).This property can be assigned to the effect of finite momentum fermion pairs scattered on the potential G k,q introduced ad hoc in (27).Physically this means that only a certain fraction of the fermion states is expelled from the vicinity of the Fermi surface.Figure 4 shows both parts of the spectral function for the pairing potential (35) with n = 100, which yields the strongest suppression of the low lying fermion states plotted in figure 3. Two modes of the coherent part have similar spectral weights as the BCS coherence factors.However, in distinction from the standard BCS solution, the corresponding quasiparticle dispersion ±ξ k (∞) in this case is gapless.This is evidently caused by coupling to the finite momentum fermion pairs.As regards the background part, it builds up mainly near the Fermi surface k F .Such fermion states are located around the quasiparticle energies ±ξ k (∞).We estimated that for k = k F , these states contribute nearly 30 percent of the total spectral weight.
The off-diagonal single particle Green's function ĉk↑ ; ĉ−p↓ ω has a structure similar to (43).Skipping unnecessary technical details we show the resulting expression for the expectation value of the order parameter We numerically investigated the momentum dependence of the order parameter ĉq−k↓ ĉk↑ and we noticed its rather fast decrease against |q|.For n = 10, 50 and 100, we obtained that the magnitude of the order parameter at |q| = 0.2π/a, which is nearly a hundred times smaller than ĉ−k↓ ĉk↑ .This indicates that finite momentum fermion pairs are less favored in the system.

Coupling to the boson mode
In this section we apply the same method to a non-trivial problem describing the itinerant fermions coupled to the dispersionless boson mode The boson mode can be regarded, for instance, as a pairing field derived from the Hubbard Stratonovich transformation in the system of correlated fermions.The same model (45) can also describe the Andreev tunneling between c-fermions and some pair reservoir denoted by b-particles.
Another possibility is to think of the case where itinerant fermions (e.g.conduction band electrons) coexist and interact with some localized fermion pairs [10].The model ( 45) is also frequently applied to the description of the ultracold fermion atoms coupled to the weakly bound molecules effectively leading to the resonant Feshbach scattering [15].
To keep a conserved total number of particles Ntot = k,σ ĉ † kσ ĉkσ + 2 b †b we apply the grand canonical ensemble.Symbol Ω stands for the boson energy measured from 2µ and g k denotes the coupling constant.It can be shown that physical properties of this model depend only on the magnitude of g k .In other words, all |g k |e iφ k lead to identical results independently of the phase φ k .
In the mean field treatment of the Hamiltonian (45), one usually introduces the linearization bĉ [10].This idea relies on the assumption that there exists a finite amount of the Bose Einstein (BE) condensation of b particles.Hamiltonian (45) is simplified to the ordinary BCS problem, where ∆ ≡ −g b .However, let us notice that the BE condensation of infinitely heavy boson (a discrete energy level Ω means that the effective mass m = ∞) is not allowed.Our present analysis based on the continuous unitary transformation shows a route to go beyond such a mean field approximation.
To construct the continuous transformation decoupling c from b particles we choose the following generating operator Applying (46) to the flow equation ( 2) we obtain the following l-dependent Hamiltonian [16] Ĥ(l) with U k,k (0) = 0 and the higher order terms involving density-density interaction between c and b objects are linearized through the normal ordering [1].In a standard way we derive the following set of flow equations where n F kσ denotes the fermion occupancy and n B = b †b refers to the statistically averaged boson number.From the formal equation (50 To have a complete set of the flow equations ( 48)-(48) we need the distribution function n F kσ .From the flow equation ( 12) we determine the following parametrizations [18] ĉk↑ where l-dependent coefficients satisfy with the boundary condition u k (0) = 1 and v k (0) = 0.In equations (54,55) we neglected the terms proportional to g k (l)U k1,k2 (l) which can eventually contribute some higher order corrections, of the order ∼ g 3 k .Combining ( 54) and ( 55) we obtain the invariance which guarantees that operators ĉ( †) kσ (l) defined in (52,53) obey the anticommutation relations.Various expectation values Ô are easy to carry out in the limit l → ∞ because fermions are decoupled from the boson field.Some delicate problem causes the two-body interaction U k,k (∞).As will be shown below this potential is small, so in the lowest order perturbation theory we incorporate its effect via the Hartree shift U k,k (∞)n F kσ to fermion energies.However, this weak point of the present analysis can be systematically improved.Figure 5 confirms that the two-body potential is indeed weak.Using the parameterization (52,53) we find the normal single particle Green's function with quasiparticle energies due to the invariance (56).Finally, the momentum distribution function is given by We numerically explored the set of coupled flow equations ( 48 2 Ω(∞).We used for computations g k = 0.05D, Ω/2 + µ = −0.3D,kBT = 0.015D.
In figure (6) we show the single particle fermion spectrum which turns out to be gapped.Two branches of the excitations E k,ν are discontinuous around the energy Ω(∞)/2 (the dashed line).We have checked that the value of the gap is rather independent of temperature.The corresponding spectral weights are shown in figure 7.They have a behavior similar to the BCS formfactors but we notice only a negligible dependence on temperature.Let us emphasize that the gapped spectrum shown in figure 6 is not related to any off-diagonal order parameter.According to the parameterization (52, 53) we determine the anomalous Green's function ĉ † −k↓ ; c † k↑ and we find it to be identically zero.The above mentioned structure should be referred to as a pseudogap.
The broken symmetry can arise if and only if a certain fraction of b particles becomes BE condensed [19].There are two possibilities for this to occur: a) either we assume bosons to be mobile right from the outset of the problem, or b) we allow for a finite momentum exchange q = 0 in the interaction term ĉ † k↑ ĉ † q−k↓ bq .The second option has been explored in the literature by several groups using various manybody techniques (for a representative list of the references see the review paper [20]).The purpose of our present study was to show that the temporal quantum fluctuations alone can induce the gapped excitation structure without involving any spatial fluctuations of the order parameter.In the real systems there would always exist both temporal and spatial fluctuations and the latter cause that upon lowering the temperature, the pseudogap smoothly evolves into the true gap of superconducting state [21].

Conclusions
By means of the flow equation method [1] we analyzed the bilinear Hamiltonians describing fermion systems having pairing interactions.This new technique becomes more and more popular [3] owing to its conceptual simplicity which allows us to go beyond the framework of various perturbative approximations [22].This method is based on the continuous canonical transformations in the course of which all the parameters are gradually renormalized.In some way this can be compared to the renormalization of the coupling constants within the RG procedure [4].The high energy excitations are renormalized here mainly in the beginning, while the low energy excitations should be dealt with until the end of transformation.
The first part of our paper illustrates the way of reproducing the rigorous solution of the BCS model using a continuous version of the Bogoliubov transformation [6].In particular, we show how one can approach the symmetry broken problems which usually is not an easy task for the RG methods [7,8].This scheme can be extended to the analysis of the two-body interactions where various kinds of instabilities can arise.Some results for the 2-dimensional Hubbard model using the flow equation method have been already reported [23].
In the second part we focused on the case where fermion pairs are coupled to some single level boson field.Such a situation can take place when the correlated fermion system is affected by bosonic modes originating from the quantum fluctuations in the systems with a reduced dimensionality like the high temperature superconductors [24].This sort of physics has recently been intensively studied for the ultracold fermion atoms where the interaction with the weakly bound molecules leads to the Feshbach resonance and turns out to be a driving mechanism for the atomic superfluidity [25].
We selfconsistently solved the corresponding set of flow equations using the numerical Runge Kutta algorithm.We showed that the interactions are responsible for the appearance of a gap in the single particle fermion spectrum.This gap is centered around the boson energy 1  2 Ω(∞) and does not signify any symmetry breaking because no order parameter is present in the system.The offdiagonal order parameter can eventually appear if the bosons have a finite mobility (characterized by some dispersion Ω q different from the discrete energy Ω considered here) when those bosons undergo the BE condensation [19,16].Here, we emphasize that a single particle gap can appear in a normal state solely due to the temporal quantum fluctuations.In the realistic systems with additional spatial fluctuations, such a normal state pseudogap can be expected to smoothly evolve into the true gap of the superconducting state.

Acknowledgement
The author kindly acknowledges F. Wegner, J. Ranninger and K.I.Wysokiński for valuable discussions.This work is partly supported by the Polish Ministry of Science and Education under the grant NN202187833.

Figure 1 .
Figure 1.Renormalization of ∆ k (l) and ξ k (l) during the flow for several values of l as denoted in the legend.At the starting point l = 0, we used the isotropic coupling ∆ k (0) = ∆.Energies are expressed in units of ∆ while l is expressed in units of ∆ −2 .

Figure 2 .
Figure 2. Flow of the coherence factor u 2 k (l) for three fermion states of the initial energy ξ k = 0.01, 0.1 and 1 (in units ∆).Inset shows the effective value of the coherence factor u 2 k (∞) versus ξ k .

Figure 3 .
Figure 3.The local density of states ρ(ω) = ¡ k A(k, ω) of fermions coupled to the potential G k,q given in (35) with n = 10, 50 and 100.The density is normalized to ρ0(ω) of the nointeracting system.

Figure 5 .
Figure 5.The resonant-like scattering potential U k,k (∞) of the interactions induced between fermions by eliminating the boson mode.

2 Figure 7 .
Figure 7.The spectral weights corresponding to the energy branches shown in figure 6.