text:
# Observation of breather-like states in a single Josephson cell

## I Introduction

Nonlinearity and discreteness in lattices lead to the existence of spatially localized stable solutions aubry ; takeno ; flach . These solutions, called discrete breathers (DBs), are characterized by time periodicity and a strong spatial localization of energy, on the scale of the lattice constant. The localized mode can be vibrational or rotational, depending on the lattice and type of excitation. In the second case the DB is called *rotobreather*. DB were found to exist in a large variety of nonlinear discrete lattices including weakly coupled optical wave guides eisenberg , crystal lattice in solids swanson , antiferromagnets schwarz  and Josephson junction arrays trias ; peter ; peter2 ; trias2 ; marcus . The latter system is a convenient experimental tool to investigate the existence and dynamics of DBs. The discrete lattice here is constituted by an array of coupled underdamped small Josephson junctions (JJs).

A JJ is composed by two superconducting layers separated by a thin tunnel barrier. Its dynamics can be described in the resistively-capacitively-shunted-junction model (RCSJ) bar . The equation for the superconducting phase difference between the two superconductor layers constituting the junction, $`\phi `$, is written in the normalized form as

$`\gamma =\ddot{\phi }+\alpha \dot{\phi }+\mathrm{sin}\phi .`$ (1)
Here the time unit is normalized by the inverse of the Josephson plasma frequency $`\omega _p`$, $`\gamma `$ is the bias current $`I`$ normalized by its critical value and $`\alpha `$ is the damping. This equation is formally equivalent to the equation of motion of a weakly damped pendulum in presence of an external torque $`\gamma `$. This mechanical analog allows for an easy description of the two possible states of the JJ. If we apply a torque $`\gamma >1`$ the pendulum will start to rotate with an average frequency $`<\dot{\phi }>`$. After the pendulum is driven to the rotating state, it will stay in this state, even for $`\gamma <1`$, as long as the external torque is large enough to compensate the damping term. So that at some range of $`\gamma `$ both static and rotating states are possible, depending on the initial conditions. By the second Josephson equation $`V=\mathrm{\Phi }_0/2\pi <\dot{\phi }>`$, these two states are the equivalent, respectively, to the superconducting and resistive states of the JJ. In experiment by measuring the voltage on the JJ it is possible to distinguish between these two states.

Josephson ladders allow for the existence of both oscillatory and rotational localized modes floria . However only the latter ones are easy to detect experimentally. A rotobreather in the Josephson ladder is constituted by some of the JJs in the resistive (rotating) state, while the rest of JJs stays in the superconducting (oscillating) state. The rotobreathers were experimentally detected by measuring the $`IV`$ characteristics of the JJs in different points of the ladder trias ; trias2  and also directly visualized by low temperature scanning laser microscopy peter ; peter2 .

Josephson ladders can support linear cavity modes that can interact with the localized modes trias2 ; marcus ; andrey ; misha2 . The interaction of the linear modes of the ladder with the localized modes was found to lead to nonlinear features in the characteristic of DBs. Their complex dynamic in the region of the linear cavity resonances was investigated experimentally by the measurements of the $`IV`$ characteristics of various DBs marcus .

Single Josephson cells have been proposed mazo ; abdel  to provide an insight into complex discrete breather states of larger systems such as Josephson ladders and two dimensional Josephson junction arrays. Broken-symmetry states with characteristics similar to those of discrete breathers were proposed to exist mazo ; abdel ; misha . We will call in the following these broken-symmetry states as *breather-like* states.

Superconducting loops with three Josephson junction have also been recently proposed for the realization of a qubit orlando ; yukon . Such a loop has an important property required for the qubit, namely two closely positioned minima in the dependence of its energy on the external magnetic field yukon2 . Note, however, that in this work we are interested into phase-rotational states, which are clearly different from oscillatory states discussed for qubits.

In this paper, we study a single Josephson cell with three small underdamped Josephson junctions. Two junctions of same area placed along the direction parallel to the uniform bias current, V1 V2, are called vertical junctions. The third junction located in one of the two transverse branches of the cell, is called horizontal H, see Fig. 1.

A parameter of cell *anisotropy*, $`\eta `$, is defined as the ratio between the critical current of the horizontal junction and the one of the vertical junctions, $`\eta =I_{CH}/I_{CV}`$. The parameter of anisotropy describes the coupling between the two vertical junctions, so that in the limit $`\eta 0`$ the two vertical junctions are entirely decoupled, while the larger is $`\eta `$ the stronger is the coupling. The cell is called isotropic if all three junctions have same parameters, anisotropic otherwise. The system is uniformly biased by two equal currents $`I`$, applied through the vertical branches, 1-2 and 3-4, Fig. 1. As long as the uniform bias current $`I`$ is smaller than the critical current $`I_{CV}`$, all the junctions stay in the superconducting state, Fig. 2a. As the current $`I`$ exceeds the critical value $`I_{CV}`$, the system switches to the homogeneous whirling state, constituted by both vertical junctions in the resistive state and the horizontal one in the superconducting state, Fig. 2b. Both these cases are symmetric states and there is no current on the horizontal junction, so that the two vertical junctions behave identically. We are interested in the last two cases, Fig. 2c and d, which are the inhomogeneous breather-like states. These two states have the same behavior and are characterized by the horizontal junction and one of the vertical junctions in the resistive state, while the other vertical junction remains in the superconducting state. In these cases, obviously, some part of the bias current flows through the transverse branch.

The procedure to generate breather-like states in experiments is similar to that used for Josephson ladders trias ; peter . We bias the system by a local current, $`I_{loc}`$, through the transverse branch 4-2, Fig. 1, and rise this current until the system switches to a breather-like state. Afterwards, we increase the bias current $`I`$ and simultaneously decrease $`I_{loc}`$ to zero, keeping the system in the generated state. Alternative procedure is first apply a homogeneous bias current $`I<I_{CV}`$ at a value where we expect to find the breather-like state stable, then rise the local current $`I_{loc}`$ to generate it, and finally reduce $`I_{loc}`$ back to zero. By using the two different directions of the initial local current it is possible to generate both breather-like states shown in Fig. 2.

Different states of the system can be experimentally detected by measuring the dc voltage on Josephson junctions. As for Josephson ladders, the broken-symmetry states are the consequence of the characteristic property of small underdamped Josephson junctions, that is the coexistence, for some ranges of the bias current, of two different stable states, superconducting and resistive.

## II Dynamics of the system

The dynamics of the system are described by the phase differences across the three Josephson junctions, $`\phi _{V_1}`$, $`\phi _{V_2}`$, $`\phi _H`$. Using the current conservation law in each node of the cell and the flux quantization for the cell, one can get the equations of motion for the system abdel ; misha :

$`𝒩(\phi _{V_1})`$ $`=`$ $`\gamma {\displaystyle \frac{1}{\beta _L}}(\phi _{V_1}\phi _{V_2}+\phi _H+2\pi f)`$
$`𝒩(\phi _{V_2})`$ $`=`$ $`\gamma +{\displaystyle \frac{1}{\beta _L}}(\phi _{V_1}\phi _{V_2}+\phi _H+2\pi f)`$ (2)
$`𝒩(\phi _H)`$ $`=`$ $`{\displaystyle \frac{1}{\beta _L\eta }}(\phi _{V_1}\phi _{V_2}+\phi _H+2\pi f),`$
where the time is normalized by the inverse of the Josephson plasma frequency $`\omega _p`$. The influence of the external magnetic field is introduced by the term $`2\pi f/\beta _L`$. The parameters of the system are the uniform *bias current* normalized by the critical current of the vertical junctions, $`\gamma =I/I_{CV}`$, the normalized *self-inductance* of the cell, $`\beta _L=2\pi LI_{CV}/\mathrm{\Phi }_0`$ (where $`L`$ is its real inductance), the *damping*, $`\alpha =\sqrt{\mathrm{\Phi }_0/(2\pi I_CR^2C)}`$, and the frustration parameter $`f=\mathrm{\Phi }_{ext}/\mathrm{\Phi }_0`$, that can be tuned by the external magnetic field $`H`$. The operator $`𝒩`$ represents the current through a single junction in the RCSJ model and is defined as: $`𝒩(\phi )=\ddot{\phi }+\alpha \dot{\phi }+\mathrm{sin}\phi `$.

Using a simplified dc model abdel , where the junctions in the resistive state are considered as resistances and the ones in the superconducting state just as shorts, one can estimate the dimensionless breather frequency $`\mathrm{\Omega }=<\dot{\phi }>`$, that is in fact the voltage across a junction in the resistive state normalized by $`\mathrm{\Phi }_0\omega _p/2\pi `$:

$`\mathrm{\Omega }={\displaystyle \frac{\gamma }{\alpha (1+\eta )}}.`$ (3)
From that one can also evaluate the region of existence of the breather-like state abdel . The minimum current at which a broken-symmetry state may exist is

$`\gamma _r={\displaystyle \frac{4\alpha }{\pi }}(1+\eta ).`$ (4)
At $`\gamma <\gamma _r`$ the system is retrapped to the superconducting state. The maximum current for a breather-like state is

$`\gamma _c={\displaystyle \frac{1+\eta }{1+2\eta }}.`$ (5)
Above this current the system switches to the homogeneous rotating state.

The simplified dc model neglects the influence of the small oscillations of the Josephson phase that may cause an instability and nonlinearity in the characteristic of the breather-like state. Neglecting the damping and linearizing the equation of motion, Eqs. (2), around the breather-like state, Eq. (3), it is possible to calculate the two characteristic frequencies of these oscillations abdel ; misha :

$`\omega _\pm =\sqrt{F\pm \sqrt{F^2G}},`$ (6)
where: $`F`$ $`={\displaystyle \frac{1}{2}}\mathrm{cos}(c_1)+{\displaystyle \frac{1+2\eta }{2\eta \beta _L}}`$
$`G`$ $`={\displaystyle \frac{1+\eta }{\eta \beta _L}}\mathrm{cos}(c_1).`$
Here $`c_1`$ is the phase shift of the vertical junction in the superconducting state evaluated by the dc model, $`c_1=\mathrm{arcsin}[\gamma (1+2\eta )/(1+\eta )]`$. The presence of the electromagnetic oscillations leads then to *primary* resonances in the breather-like state, when the frequency of the breather $`\mathrm{\Omega }`$ locks to the frequency of the electromagnetic oscillations $`\omega _+`$ or $`\omega _{}`$, and to *parametric* and *combination* resonances, respectively, at $`\mathrm{\Omega }=n\omega _\pm `$, where $`n2`$ is an integer, and $`\mathrm{\Omega }=n(\omega _+\pm \omega _{})`$, where $`n`$ is an integer abdel ; misha .

Using the rotational approximation kulik , two different dependencies of the resonant step amplitude $`\mathrm{\Delta }\gamma `$ on the frustration have been obtained for *primary* and *parametric* resonances misha :

$`\mathrm{\Delta }\gamma _{\omega _\pm }`$ $``$ $`\left|\mathrm{sin}\left({\displaystyle \frac{c_1}{2}}+{\displaystyle \frac{\beta _L\gamma \eta }{2(1+\eta )}}+\pi f\right)\right|`$ (7)
$`\mathrm{\Delta }\gamma _{n\omega _\pm }`$ $``$ $`\sqrt{(1\eta )^2+4\eta \mathrm{cos}^2\left({\displaystyle \frac{c_1}{2}}+{\displaystyle \frac{\beta _L\gamma \eta }{2(1+\eta )}}+\pi f\right)}.`$
Expressions (7) and (LABEL:para) have the same periodicity in $`f`$ as the critical current in the whirling state, but, depending on the parameters of the system, they may be asymmetric with respect to $`f=0`$. Moreover, the external magnetic field can considerably change the amplitude of the resonance and, in the case of primary resonance, even make it to vanish for some values of the frustration.

## III Experimental results

Our experimental system is a single cell with the hole area $`4\times 4`$ $`\mu `$m2, containing three small underdamped Nb/Al-AlOx/Nb Josephson junctions hypres , see Fig. 1. The critical current density, measured at 4.2 K, is $`J_C=133`$ A/cm2. The area of each vertical junction is about $`7\times 7`$ $`\mu `$m2, while the horizontal one is $`3\times 4`$ $`\mu `$m2. The anisotropy of the system, evaluated from the design, is then $`\eta =0.24`$. The parameter of self-inductance of the cell at $`4.2`$ K, evaluated by measurement of a SQUID of same geometry, is $`\beta _L=1.47`$. The two equal bias currents are introduced in parallel using two large external resistances of $`1.5`$ k$`\mathrm{\Omega }`$, so that the bias current remains uniform also in the broken-symmetry states. The damping of the system at 4.2 K is evaluated from the retrapping current of the homogeneous whirling state, $`\alpha (\pi /4)(I_r/2I_{CV})`$, that yields $`\alpha 0.03`$. The measurements are performed in a 4He cryostat cryo  in the temperature range between 4 K and 8 K, so that the parameters of the system $`J_C`$, $`\alpha `$, $`\beta _L`$ and $`\omega _p`$ can be varied. For each temperature the damping $`\alpha `$ and the critical current density $`J_C`$ are evaluated from the characteristic of the homogeneous whirling state, while the parameter of self inductance $`\beta _L`$ is evaluated from the critical current density. Finally, the value of the plasma frequency is obtained at every temperature from fitting the voltage of the parametric resonant step to the expected value $`2\omega _+`$.

For convenience of comparing with theory, the experimental data are shown in the previously introduced normalized units of the bias current $`\gamma =I/I_{CV}`$, breather frequency $`\mathrm{\Omega }=V/(\omega _p\mathrm{\Phi }_0/2\pi )`$ and frustration $`f=\mathrm{\Phi }_{ext}/\mathrm{\Phi }_0`$.

In the range of temperatures from 4.2 K to 7.75 K ($`0.02\alpha 0.35`$) it was possible to excite both breather-like states and they were found stable in a finite range of the bias current. No breather-like state was found stable for temperatures above 7.75 K ($`\alpha >0.35`$).

The minimum bias current below which the breather is retrapped increases linearly with the damping, see Fig. 3, in good agreement with the theoretical expectation of Eq (4). The maximum current that allows for the existence of the breather-like state slightly decreases with increasing the damping, see Fig. 3, whereas from Eq. (5) it is expected to be constant ($`\gamma _c=0.84`$).

In fact, the temperature variation leads to an alteration of several parameters, see Fig. 4. In particular the decrease of the self inductance $`\beta _L`$ at high temperatures leads to stronger interactions of the breather-like state with the small oscillations. This effect may strongly affect (in the range $`0.1<\beta _L<1`$) the region of existence of the breather-like, especially its upper border abdel . In fact, the experimental data show deviation from the expected value for $`\alpha 1.15`$, which is in the range of the self inductance between 0.9 and 0.3. No resonances were observed in the $`IV`$ characteristics at temperatures below 6 K. The small damping of the junctions at this temperatures seems to be unfavorable for reaching the small frequencies of the electromagnetic oscillations where the resonances are expected ($`\mathrm{\Omega }_{min}5.5`$). Increasing the temperature of the system had two helpful effects for seeing the resonances, namely the increase of damping and decrease of self inductance, see Fig. 4. As the consequence of the larger damping, the breather-like states were stable down to smaller frequencies, whereas a smaller inductance moved the frequencies of the resonances up.

At the temperature of $`T=6`$ K it was possible to observe the first resonant step. The parameters of the system at this temperature are $`\alpha 0.11`$, $`J_c=92`$ A/cm2, $`\beta _L=1.02`$ and $`\omega _p/2\pi =37`$ GHz. This resonance is a step in the characteristic, shown in Fig. 5. Comparing this resonance with that at higher temperatures and with the theoretical expectation we deduce it to be the parametric resonance at $`2\omega _+`$.

At the temperature of $`T=6.65`$ K we detect two additional resonances in the characteristic of the breather-like state, see Fig. 6. At this temperature the parameters of the system are $`\alpha 0.18`$, $`J_c=70`$ A/cm2, $`\beta _L=0.77`$ and $`\omega _p/2\pi =30`$ GHz. The two larger resonances, A and B, are different by the integer factor of two, as expected for $`2\omega _+`$ and $`\omega _+`$.

The dependencies of the current amplitude on the external magnetic field for steps A and B are strongly different. The resonance A has its maximum amplitude close to zero frustration, see Fig. 6a, and disappears in a wide range around $`f=0.5`$, see Fig. 7. An other resonance B is always present but has nearly opposite behavior. At the top of the resonance B, around frustration $`f=0.4`$ or $`f=0.6`$, the system shows an instability by switching either to the homogeneous whirling state or to the superconducting state, as shows Fig. 7. In this range of frustration the measurements of the resonant step B were performed with less accuracy and faster because of its poor stability. We argue the larger step B is due to the primary resonance at $`\mathrm{\Omega }=\omega _+`$, whereas the step A is due to the parametric resonance at $`\mathrm{\Omega }=2\omega _+`$. The third smaller nonlinearity, C, is at a frequency between the other two, close to the theoretically expected value for the composed resonance at $`\mathrm{\Omega }=\omega _++\omega _{}`$, see Fig 6. This resonance was only present in a small range of frustration around zero, where from the theory, the composed resonance is expected to be maximum. Because of its reduced amplitude this resonant step has no hysteresis and its dependence on the magnetic field was not measured.

In Fig. 8b, c we present measurements of the current amplitude of the resonant steps A and B on the external magnetic field. The parametric step at $`T=6`$ K, see Fig. 5, is higher than that at $`T=6.65`$ K, Fig. 6. Moreover, because of the larger damping in the latter case the step vanishes in a larger range of the frustration ($`0.95<f<0.25`$ and $`0.05<f<0.75`$). Therefore, in order to observe a more clear modulation of this resonance, we report in Fig. 8c its dependence on the external magnetic field at $`T=6`$ K.

As expected from the theory, Eq. (7) and Eq. (LABEL:para), the magnetic field dependence of the resonance amplitudes show an asymmetry with respect to $`f=0`$. Note, that such asymmetry is not observed for $`I_C(f)`$, see Fig. 8(a). The asymmetry depends on the parameter of the system and particularly on the bias current where the resonance occurs. This behavior can be qualitatively attributed to the current flowing in the transverse branch of the cell in the broken-symmetry state. The experimental data are in rather good agreement with the theoretical expectations. We fitted the maximum and minimum of the resonant steps in order to find the proportionality coefficient in Eqs. (7) and (LABEL:para). The small variations of the experimental data around the maximum of the primary resonant step could be due to the appearance of instabilities, as discussed in relation to Fig. 7. This behavior was also found in simulations misha . In the case of parametric resonance the presence of a range of frustration, that increases with the damping $`\alpha `$, where resonant step disappears, is also expected misha . Note, that the region where the experimental data show the largest deviation from the expected behavior is nearly the same for the two cases ($`0.3<f<0.6`$).

## IV Conclusions

We observed breather-like states of broken symmetry in an anisotropic single-cell system. For temperatures lower than 6 K the current-voltage characteristics of the breather-like states were stable and did not show any resonances. Increasing the dissipation by rising the temperature of the system we decreased the frequencies of rotation of the breather-like state in such a way that it matched the frequencies of the resonances in the cell. In this region, the primary, composed and parametric resonances were observed. As expected, in all three cases, the interaction of the electromagnetic oscillations with the breather-like state leads to resonant steps in the current-voltage characteristic. The large steps due to the primary and parametric resonances show hysteresis, whereas the composed one induces just a small nonlinearity in the characteristic. The measured dependence on the external magnetic field of the magnitude of the primary and parametric resonant steps is in good agreement with theory. In the case of primary resonance, the external frustration can strongly increase the amplitude of the resonant step and lead to instability of the breather-like state. The composed resonance, for this values of the damping, was present just in a small range around zero frustration and did not lead to instabilities. The region of existence of the breather-like state was also measured and found in good accordance with the expectations. At smaller values of the self inductance parameter, the stronger interaction of the breather with the electromagnetic oscillations leads to a contraction of the region of existence.

## V acknowledgments

We would like to thank A. Benabdallah, P. Binder, M. V. Fistul, S. Flach, A. E. Mirishnichenko and M. Schuster for the fruitful discussions and for sending us their latest works before publication. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) and by the European Union under RTN project LOCNET HPRN-CT-1999-00163.


format: text

text:
# Spin collective mode and quasiparticle contributions to STM spectra of 𝑑-wave superconductors with pinning

## 1 Introduction

The observation of a halo of checkerboard modulation in the local density of states (LDOS) around magnetic field-induced vortex cores by Hoffman et al.  in optimally doped Bi2Sr2CaCu2O8+δ has focused attention on the nature of the spin and charge collective modes in the cuprate superconductors. Charge order in this region was anticipated theoretically . By ‘charge order’ we mean here periodic spatial modulations of all local observables invariant under spin rotations and time-reversal, such as the exchange and pairing energies, and the LDOS itself. If the modulation is primarily on quantities defined on a bond, rather than in the charge density on a site, then such order is better denoted as ‘bond order’; it was this bond order which was theoretically anticipated . The general connection relating such charge/bond order to the ‘square’ of collinear, non-bipartite, spin correlations  (see Section 2 below) was used to develop a model  relating the spatial structure and energy spectrum of the charge order to the dynamic spin fluctuation spectrum as measured by neutron scattering . There have also been other discussions of competing order within vortex cores, including enhanced dynamic antiferromagnetism in a spin-gap state , static Néel order , and other models with static spin density wave order .

More recent measurements of Howald et al.  have shown that the LDOS modulations in the vortex halo also persists in zero magnetic field, albeit with a much weaker amplitude: here the pinning of the charge/bond order is presumably due to impurities and other imperfections in the Bi2Sr2CaCu2O8+δ crystal. Theoretical modeling of the STM spectrum has indicated  that the order is likely to be primarily on bond variables, as may be derived from a spin-Peierls state , i.e. a bond order wave . This is in contrast to the prediction of a charge density modulation on the Cu sites in the anti-phase domain wall picture of ‘stripes’ , which is not observed .

Hoffman et al.  have also reported related measurements in zero magnetic field, and have focused attention on modulations of the LDOS at wavevectors throughout the Brillouin zone, and on their dispersion as a function of energy. They have identified some wavevectors with quasiparticle excitations of the superconductor scattered by random impurities, as suggested in Refs..

The purpose of the present paper is to discuss the relationship between the quasiparticle and collective mode contributions to the LDOS modulations in a $`d`$-wave superconductor. We will use precisely the model proposed earlier by Polkovnikov et al.  for the pinning of the spin collective mode: we show that the earlier expressions for the induced LDOS modulations also contain additional contributions which can be interpreted as those due to quasiparticle scattering, with a spectrum related to observations . The relative strength of the collective and quasiparticle contributions can be varied by changing the gap to the spin collective mode: experimentally the gap can be reduced (and the collective mode contribution enhanced) by lowering the doping or by applying a magnetic field , and this should allow a test of our predictions. Some related observations have also been made independently by Han . Our main conclusion is that the LDOS modulations observed by STM contain contributions from both the spin collective mode and the quasiparticles: as described in our simple model , these are relatively independent low energy excitations of the BCS superconductor, but they can become strongly coupled in some sectors of the Brillouin zone if permitted by energy and wavevector conservation ; they should not be viewed as mutually exclusive sources of LDOS modulations.

## 2 Model

We assume the cuprate superconductor has a dominant spin collective mode near the wavevectors $`𝐊_x=(3\pi /4,\pi )`$ and $`𝐊_y=(\pi ,3\pi /4)`$, and so write for the spin operator on the lattice site $`𝐫`$ at imaginary time $`\tau `$

$$S_\alpha (𝐫,\tau )=\text{Re}\left[e^{i𝐊_x𝐫}\mathrm{\Phi }_{x\alpha }(𝐫,\tau )+e^{i𝐊_y𝐫}\mathrm{\Phi }_{y\alpha }(𝐫,\tau )\right],$$
(1)
where $`\mathrm{\Phi }_{x,y\alpha }`$ are the spin density wave (SDW) order parameters which are assumed to be smooth functions of spacetime. As an aside, and for completeness, we mention that charge/bond order parameters can also be constructed out of the $`\mathrm{\Phi }_{x,y\alpha }`$. We define

$`Q_𝐚(𝐫,\tau )`$ $``$ $`S_\alpha (𝐫,\tau )S_\alpha (𝐫+𝐚,\tau )`$ (2)
$``$ $`\text{Re}\left[e^{2i𝐊_x𝐫+i𝐊_x𝐚}\mathrm{\Phi }_{x\alpha }^2(𝐫,\tau )\right]+\mathrm{}`$
where $`𝐚`$ is a vector representing a bond (say the nearest neighbor vector), we have assumed that $`\mathrm{\Phi }_{x\alpha }`$ does not vary significantly over the spatial distance $`𝐚`$, and the ellipses denote numerous other similar terms which can be deduced from (1). For $`𝐚`$ a nearest-neighbor vector, $`Q_𝐚`$ measures the modulations in the exchange energy, while for $`𝐚=0`$, $`Q_𝐚`$ measures the local charge density (for the $`t`$-$`J`$ model the charge density is linearly related to the on-site $`S_\alpha ^2`$). Thus we see from (2) that, as noted earlier, the charge/bond order parameter is proportional to the square of the SDW order parameter. In a similar manner, we will see below that the LDOS modulations are also sensitive to the square of the SDW order and so occur at wavevectors $`2𝐊_{x,y}`$: however our computation will include retardation effects missing from a phenomenological correspondence analogous to (2). All our computations below will be in a regime where

$$S_\alpha =\mathrm{\Phi }_{x\alpha }=\mathrm{\Phi }_{y\alpha }=0$$
(3)
(because spin rotation invariance is preserved), but we will nevertheless have the rotationally invariant

$$\underset{\alpha }{}\mathrm{\Phi }_{x,y\alpha }^20,$$
(4)
because translational invariance is broken by a pinning potential. This should be contrasted with the distinct physical pictures of Refs.  all of which break spin rotation invariance with $`S_\alpha 0`$.

To keep this paper self-contained, we now briefly describe the model used to compute the LDOS modulations; for details see Ref. . The effective action for the $`\mathrm{\Phi }_{x,y\alpha }`$ fluctuations is assumed to be Gaussian:

$$𝒮_\mathrm{\Phi }=T\underset{𝐪,\omega _n,\alpha }{}\chi _x^1(𝐪,\omega _n)|\mathrm{\Phi }_{x\alpha }(𝐪,\omega _n)|^2+(xy)$$
(5)
where $`\omega _n`$ are Matsubara frequencies. The $`\chi _{x,y}(𝐪,\omega )`$ denote the dynamic spin susceptibility as a function of wavevector $`𝐪`$ (measured from the ordering wavevectors $`𝐊_{x,y}`$) and frequency $`\omega `$, which is measured by inelastic neutron scattering. In the absence of available experimental input for Bi2Sr2CaCu2O8+δ, we will continue to use the simple phenomenological form  expected near a SDW ordering quantum phase transition in a superconductor: $`\chi ^1(𝐪,\omega _n)=\omega _n^2+c^2𝐪^2+\mathrm{\Delta }^2`$, where $`c`$ is a spin wave velocity, and $`\mathrm{\Delta }`$ is the spin gap to collective excitations near the wavevectors $`𝐊_{x,y}`$ (we have assumed $`\chi _x=\chi _y\chi `$ for simplicity). Han  has used the same model, but with microscopic computations on the lattice $`t`$-$`J`$ model as input for $`\chi `$.

The all-important ingredient which produces an observable static modulation in the STM signal is the pinning of the sliding degree of freedom of the SDW fluctuations. This pinning is realized by any imperfection (impurity, vortex core, etc.) which breaks the lattice translational symmetry. We employ the simplest possible pinning term which is invariant under the SU(2) spin symmetry,

$$𝒮_{\mathrm{pin}}=\underset{\alpha }{}d\tau [\zeta _x\mathrm{\Phi }_{x\alpha }^2(𝐫_0,\tau )+(xy)+c.c.],$$
(6)
where $`\zeta _{x,y}`$ are complex coupling constants representing the pinning potential.

The electronic quasiparticles of the $`d`$-wave superconductor are described by a standard BCS model:

$$H_{\mathrm{BCS}}=\underset{𝐤}{}\mathrm{\Psi }_𝐤^{}\left[(\epsilon _𝐤\mu )\tau ^z+\mathrm{\Delta }_𝐤\tau ^x\right]\mathrm{\Psi }_𝐤.$$
(7)
Here $`\mathrm{\Psi }_𝐤=(c_𝐤,c_𝐤^{})`$ is a Nambu spinor at momentum $`𝐤=(k_x,k_y)`$, $`\tau ^{x,y,z}`$ are Pauli matrices in particle-hole space, and $`\mu `$ is a chemical potential. For the kinetic energy, $`\epsilon _𝐤`$ we have first ($`t`$) and second ($`t^{}`$) neighbor hopping, and the gap function is $`\mathrm{\Delta }_𝐤=(\mathrm{\Delta }_0/2)(\mathrm{cos}k_x\mathrm{cos}k_y)`$. The full Hamiltonian for the conduction electrons includes a coupling to the collective SDW fluctuations:

$`H`$ $`=`$ $`H_{BCS}+\gamma {\displaystyle \underset{𝐫}{}}c_\mu ^{}(𝐫)\sigma _{\mu \nu }^\alpha c_\nu (𝐫)\times `$
$`\times `$ $`\left(\mathrm{\Phi }_{x\alpha }(𝐫)e^{i𝐊_x𝐫}+\mathrm{\Phi }_{x\alpha }^{}(𝐫)e^{i𝐊_x𝐫}+(xy)\right),`$
where $`\sigma ^\alpha `$ are the Pauli spin matrices, and $`\gamma `$ describes the scattering of the Bogoliubov quasiparticles of the superconductor off the SDW fluctuations.

The predictions of this model  for the quasiparticle LDOS, to second order in the quasiparticle-collective mode coupling $`\gamma `$ and to first order in the pinning strength $`\zeta `$, can be written as

$`\delta \rho (𝐫,\omega )\mathrm{Im}\mathrm{Tr}{\displaystyle \frac{1+\tau _z}{2}}\times `$
$`{\displaystyle \underset{𝐫_1,𝐫_2}{}}G_0(𝐫𝐫_1,\omega )\mathrm{\Sigma }(𝐫_1,𝐫_2,\omega )G_0(𝐫_2𝐫,\omega ),`$ (9)
where $`G_0`$ is the electronic Green’s function in the canonical translationally invariant BCS superconductor described by (7). There is a contribution to the self energy $`\mathrm{\Sigma }=\mathrm{\Sigma }_{\mathrm{sp}}`$ from the coupling to the spin collective mode, and we focus only on the part which breaks translational invariance because of the pinning of the collective mode

$`\mathrm{\Sigma }_{\mathrm{sp}}(𝐫_1,𝐫_2,\omega )12\zeta \gamma ^2{\displaystyle \frac{d\nu }{2\pi }\chi (𝐫_2,\nu )\chi (𝐫_1,\nu )}`$
$`\times G_0(𝐫_2𝐫_1,\omega \nu )\mathrm{cos}(𝐊_x(𝐫_1+𝐫_2)\delta )`$
$`+(xy).`$ (10)
The parameter $`\delta `$, which is the phase of the coupling $`\zeta `$, distinguishes site centered SDWs ($`\delta =0`$) and bond centered SDWs ($`\delta =\pi /2`$). Following the STM experiments , we will compute the Fourier transform of the LDOS:

$`\delta \rho _{\mathrm{sp}}(𝐪,\omega )=\mathrm{Tr}{\displaystyle \frac{1+\tau _z}{4}}\{\text{Im}{\displaystyle }{\displaystyle \frac{d^2k}{(2\pi )^2}}G_0(𝐤,\omega )\times `$
$`[\mathrm{\Sigma }_{\mathrm{sp}}(𝐤,𝐤+𝐪,\omega )+\mathrm{\Sigma }_{\mathrm{sp}}(𝐤,𝐤𝐪,\omega )]G_0(𝐤+𝐪,\omega )`$
$`i\text{Re}{\displaystyle }{\displaystyle \frac{d^2k}{(2\pi )^2}}G_0(𝐤,\omega )[\mathrm{\Sigma }_{\mathrm{sp}}(𝐤,𝐤+𝐪,\omega )`$
$`\mathrm{\Sigma }_{\mathrm{sp}}(𝐤,𝐤𝐪,\omega )]G_0(𝐤+𝐪,\omega )\}.`$ (11)
Note that the imaginary part of this expression vanishes in the presence of inversion symmetry, $`\delta \rho (𝐫,\omega )=\delta \rho (𝐫,\omega )`$. The results (10,11) are for a single pinning center, and the realistic many impurity case will be similar provided interference between impurities is small.

It is not difficult to see from (9-11) that in the limit of a large spin gap ($`\mathrm{\Delta }\mathrm{}`$), the modulations in the LDOS predicted by (11) are identical to those predicted in a model of quasiparticle scattering off isolated, point-like static impurities or other imperfections considered in Ref. . For a general scattering potential $`u(𝐪)`$, and impurity-induced variation in the superconducting gap $`v(𝐪_1,𝐪_2)`$, the expression for the LDOS modulation due to static impurities in the Born approximation, $`\delta \rho _{\mathrm{imp}}(𝐪,\omega )`$, is formally identical to Eq. (11) with $`\mathrm{\Sigma }_{\mathrm{sp}}`$ replaced by

$$\mathrm{\Sigma }_{\mathrm{imp}}(𝐤,𝐤+𝐪,\omega )=u(𝐪)\tau _z+v(𝐤,𝐤+𝐪)\tau _x.$$
(12)
The following section will compare the predictions of (12) to those in (11) of a pinned dynamic SDW with a finite $`\mathrm{\Delta }`$.

It is interesting to note from (11,12) that for any real scattering potential, with $`u^{}(𝐪)=u(𝐪)`$, the LDOS modulation at wavevector $`𝐪`$ is determined by $`d^2k\text{Im}\left[G_0(𝐤)G_0(𝐤+𝐪,\omega )\right]`$. Contributions proportional to $`d^2k\text{Im}\left[G_0(𝐤)\right]\text{Im}\left[G_0(𝐤+𝐪,\omega )\right]`$, which are used by Ref.  to explain the STM observations of Refs. , can only come from lifetime effects which produce a complex self energy.

## 3 Results

Sample results from the dynamic SDW model (11) are shown in Figs 1-3.

For small spin gap $`\mathrm{\Delta }`$, the LDOS shows the expected  modulations near $`𝐪=2𝐊_{x,y}`$. The magnitude of this modulation depends strongly on the value of $`\mathrm{\Delta }`$, and the strong expected dependence of $`\mathrm{\Delta }`$ on doping and applied magnetic field  ($`\mathrm{\Delta }`$ decreases with decreasing $`\delta `$ and increasing field) will lead to a corresponding variation in the LDOS peaks near $`𝐪=2𝐊_{x,y}`$. This peak also disperses as a function of $`\omega `$, with the dispersion being stronger for larger values of $`\mathrm{\Delta }`$ as shown in Fig. 3; the decrease in the $`𝐪`$ value of the peak with increasing $`\omega `$ is consistent with observations .

As is clear from Figs. 1 and 2, the LDOS also has peaks in its modulation at a number of other wavevectors. These arise from the two flanking quasiparticle Green’s functions in (11), and are consequently sensitive to all details of quasiparticle spectrum. These additional features are also present in the static impurity model (12) whose results are shown in Fig. 4

. It is instructive to compare the static impurity scattering prediction in Fig 4 to the dynamic spin collective mode predictions in Figs. 1 and 2: apart from the strong contributions at $`𝐪=\pm 2𝐊_{x,y}`$ present in the latter, the former has the same general features as the large $`\mathrm{\Delta }`$ cases of the latter. This is one of the main points of this paper.

The physics of the structure in the static impurity scattering LDOS modulation, $`\delta \rho _{\mathrm{imp}}(𝐪,\omega )`$ in (12), has already been discussed by Wang and Lee . Assuming spatial inversion symmetry, the expression (12) can be rewritten as:

$`\delta \rho _{\mathrm{imp}}(𝐪,\omega )=2{\displaystyle }{\displaystyle \frac{d^2k}{(2\pi )^2}}\delta (\omega ^2+\epsilon _𝐤^2+\mathrm{\Delta }_𝐤^2)\times `$
$`(u(𝐪){\displaystyle \frac{(\omega \epsilon _𝐤)(\omega \epsilon _{𝐤+𝐪})\mathrm{\Delta }_𝐤\mathrm{\Delta }_{𝐤+𝐪}}{\omega ^2+\epsilon _{𝐤+𝐪}^2+\mathrm{\Delta }_{𝐤+𝐪}^2}}`$
$`v(𝐤,𝐤+𝐪){\displaystyle \frac{(\omega \epsilon _𝐤)\mathrm{\Delta }_{𝐤+𝐪}+\mathrm{\Delta }_𝐤(\omega \epsilon _{𝐤+𝐪})}{\omega ^2+\epsilon _{𝐤+𝐪}^2+\mathrm{\Delta }_{𝐤+𝐪}^2}}).`$ (13)
This expression is not a convolution of the product of the density of states of quasiparticles at energy $`\omega `$ with momenta differing by $`𝐪`$, but rather the product of the density of states of one quasiparticle with the real part of the Green function of the other. It contains weak square-root divergences at special nesting wavevectors , where constant energy contours with momenta $`𝐤`$ and $`𝐤+𝐪`$ ($`\omega ^2=\epsilon _𝐤^2+\mathrm{\Delta }_𝐤^2=\epsilon _{𝐤+𝐪}^2+\mathrm{\Delta }_{𝐤+𝐪}^2`$) are tangent to each other at one of their the points of contact. Such wavevectors $`𝐪`$ form different lines in the Brillouin zone, but these are barely visible in our plots. It has been claimed  that the structure at related wavevectors could be enhanced for strong impurity scattering in a description beyond the Born approximation. Instead, in Fig 4, the main contributions for $`𝐪`$ along the Brillouin zone diagonals come from values of $`𝐤`$ and $`𝐤+𝐪`$ near the nodal points and the underlying Fermi surface near $`\pm (\pi ,0),\pm (0,\pi )`$ respectively. We emphasize that both these features of quasiparticle-scattering induced LDOS modulations also appear in our model  of a pinned dynamic spin collective mode. In particular, as we noted earlier, if $`\mathrm{\Delta }`$ is large, then the pinned SDW model is equivalent to an effective short-range static quasiparticle potential, along with additional weak modulations at $`\pm 2𝐊_{x,y}`$ (compare Figs 2, and 4). For smaller $`\mathrm{\Delta }`$, the effective quasiparticle potential created by the SDW fluctuations is dynamic, but the qualitative picture at wavevectors different from $`\pm 2𝐊_{x,y}`$ remains the same.

In conclusion, we have shown that the model of Ref  displays many of the features observed in recent STM experiments. Besides the modulations of the LDOS at wavevectors $`\pm 2𝐊_{x,y}`$ arising from the spin collective mode, we find other features related to the spectrum of the fermionic Bogoliubov quasiparticles (similar physics is also expected from a possible charge/bond order collective mode). The positions of the maxima at $`\pm 2𝐊_{x,y}`$ have a relatively small (but not negligible) dependence on the tunnelling energy $`\omega `$, while other features are more strongly dependent on $`\omega `$ and the details of the quasiparticle spectrum . We have also made predictions on the doping and magnetic field dependence of the relative contributions of the quasiparticle and spin-collective mode induced modulations.

Note added: Some related points and observations appear in .

We thank Seamus Davis, Jung Hoon Han, Aharon Kapitulnik, Steve Kivelson, Patrick Lee, and Vadim Oganesyan for useful discussions. A.P. and S.S. were supported by US NSF Grant DMR 0098226. M.V. was supported by the DFG through SFB 484.


text:
# Renormalized perturbation theory for Fermi systems: Fermi surface deformation and superconductivity in the two-dimensional Hubbard model

## I Introduction

Unrenormalized perturbation expansions of interacting electron systems around the non-interacting part of the Hamiltonian are generally plagued by infrared divergencies. Some of the divergencies are simply due to to shifts of the Fermi surface, while others signal instabilities of the normal Fermi liquid towards qualitatively different states, such as superconducting or other ordered phases. This problem is often treated by self-consistent resummations of Feynman diagrams, where a finite or infinite subset of skeleton diagrams, with the interacting propagator $`G`$ on internal lines, is summed. Symmetry breaking can be built into the structure of $`G`$ as an ansatz, and the size of the corresponding order parameter is determined self-consistently. This standard approach has been very useful in many cases. However, resummation schemes beyond first order (Hartree-Fock) require extensive numerics, since the full self-energy has to be determined self-consistently, and delicate low-energy structures cannot always be resolved. A more serious problem is the fact that self-energy and vertex corrections are not treated on equal footing in most feasible resummation schemes. This often leads to unphysical results.

In this work we will describe and apply an alternative procedure, which has been formulated already long ago by Nozières, and more recently been discussed in the mathematical literature as a way of carrying out well-defined perturbation expansions for weakly interacting Fermi systems. The basic idea is to choose an improved starting point for the perturbation expansion, by adding a suitable counterterm to the non-interacting part of the Hamiltonian, and subtracting it from the interaction part. The counterterm is quadratic in the Fermi operators and has to be determined from a self-consistency condition. In Sec. II we will describe how Fermi surface deformations and superconductivity can be treated by this method. Explicit expressions up to second order in the interaction are derived for the case of the Hubbard model in Sec. III. Results obtained from a numerical solution of the self-consistency equations in two dimensions will follow in Sec. IV. For the repulsive Hubbard model we have obtained superconducting solutions with d-wave symmetry in agreement with widespread expectations, and with recent renormalization group calculations which conclusively established d-wave superconductivity at weak coupling. In addition, for Fermi levels close to the van Hove singularity, deformations which break the square lattice symmetry occur. This confirms the recently proposed possibility of symmetry-breaking Fermi surface deformations (”Pomeranchuk instabilities”).

## II Renormalized perturbation expansion

We consider a system of interacting spin-$`\frac{1}{2}`$ fermions with a Hamiltonian $`H=H_0+H_I`$, where the non-interacting part

$$H_0=\underset{𝐤,\sigma }{}\xi _𝐤n_{𝐤\sigma }$$
(1)
with $`\xi _𝐤=ϵ_𝐤\mu `$ contains the kinetic energy and the chemical potential, while $`H_I`$ is a fermion-fermion interaction term. We are particularly interested in lattice systems, for which the dispersion relation $`ϵ_𝐤`$ is not isotropic. We consider only ground state properties, that is the temperature is zero throughout the whole article.

The bare propagator in a standard many-body perturbation expansion  around $`H_0`$ is given by

$$G_0(k)=\frac{1}{i\omega \xi _𝐤},$$
(2)
where $`\omega `$ is the Matsubara frequency and $`k=(\omega ,𝐤)`$. This progagator diverges for $`\omega 0`$ and $`𝐤𝐤_F`$, for any Fermi momentum $`𝐤_F`$, since $`\xi _{𝐤_F}=0`$. As a consequence, many Feynman diagrams diverge. A well-known singularity is the (usually) logarithmic divergency of the 1-loop particle-particle contribution to the two-particle vertex in the Cooper channel, which leads to a $`(\mathrm{log})^n`$ divergency of the n-loop particle-particle ladder diagram. This signals a possible Cooper instability towards superconductivity. Much stronger divergencies occur in diagrams with multiple self-energy insertions on the same internal propagator line, leading to non-integrable powers of $`G_0(k)`$. These singularities are due to Fermi surface shifts generated by the interaction term in the Hamiltonian.

The divergency problems and the superconducting instability can be treated by splitting the Hamiltonian in a different way, namely as

$$H=\stackrel{~}{H}_0+\stackrel{~}{H}_I,$$
(3)
where $`\stackrel{~}{H}_0=H_0+\delta H_0`$ and $`\stackrel{~}{H}_I=H_I\delta H_0`$, and expanding around $`\stackrel{~}{H}_0`$. The counterterm $`\delta H_0`$ must be quadratic in the creation and annihilation operators to allow for a straightforward perturbation expansion based on Wick’s theorem. It is possible to chose $`\delta H_0`$ such that $`\stackrel{~}{H}_I`$ does not shift the Fermi surface corresponding to $`\stackrel{~}{H}_0`$ any more, and divergencies due to self-energy insertions are removed. In the superconducting state spontaneous symmetry breaking can be included already in $`\delta H_0`$, with an order parameter $`\mathrm{\Delta }_𝐤`$ whose value on the Fermi surface is not shifted by $`\stackrel{~}{H}_I`$. We will now describe this procedure in more detail.

### A Normal state

A counterterm $`\delta H_0=_{𝐤,\sigma }\delta \xi _𝐤n_{𝐤\sigma }`$ leads to a renormalized dispersion relation $`\stackrel{~}{\xi }_𝐤=\xi _𝐤+\delta \xi _𝐤`$ in the unperturbed part of the Hamiltonian,

$$\stackrel{~}{H}_0=\underset{𝐤,\sigma }{}\stackrel{~}{\xi }_𝐤n_{𝐤\sigma },$$
(4)
and correspondingly to a new bare propagator

$$\stackrel{~}{G}_0(k)=\frac{1}{i\omega \stackrel{~}{\xi }_𝐤}.$$
(5)
The Fermi surface $`\stackrel{~}{}`$ associated with $`\stackrel{~}{H}_0`$ is given by the momenta $`\stackrel{~}{𝐤}_F`$ satisfying the equation $`\stackrel{~}{\xi }_𝐤=0`$. The Fermi surface of the interacting system is given by the solutions of the equation $`G^1(0,𝐤)=0`$. This surface coincides with the unperturbed one, corresponding to $`\stackrel{~}{H}_0`$, if the renormalized self-energy $`\stackrel{~}{\mathrm{\Sigma }}=\stackrel{~}{G}_0^1G^1`$ vanishes on $`\stackrel{~}{}`$, that is if

$$\stackrel{~}{\mathrm{\Sigma }}(0,𝐤)=0\text{for}𝐤\stackrel{~}{}.$$
(6)
This imposes a self-consistency condition on the counterterms which can be solved iteratively. For isotropic systems the shift of $`\xi _𝐤`$ can be chosen as a momentum independent constant, which may be interpreted as a shift of the chemical potential. For anisotropic systems, however, one generally has to adjust the whole shape of the Fermi surface. That this procedure really works at each order of the perturbation expansion has been shown rigorously for a large class of systems.

The shift function $`\delta \xi _𝐤`$ is uniquely determined by the self-consistency condition only on the (interacting) Fermi surface $`\stackrel{~}{}`$. For momenta away from the Fermi surface, $`\delta \xi _𝐤`$ can be chosen to be any sufficiently smooth function of $`𝐤`$ which does not lead to artificial additional zeros of $`\stackrel{~}{\xi }_𝐤`$.

The perturbation expansion of the renormalized self-energy $`\stackrel{~}{\mathrm{\Sigma }}`$ involves two types of vertices: the usual two-particle vertex given by the interaction $`H_I`$ and a one-particle vertex due to the counterterm $`\delta H_0`$ in $`\stackrel{~}{H}_I`$. In Fig. 1 we show the Feynman diagrams contributing to $`\stackrel{~}{\mathrm{\Sigma }}`$ up to second order in the interaction.

The above-mentioned divergencies of Feynman diagrams with self-energy insertions on internal propagator lines are removed in the renormalized expansion around $`\stackrel{~}{H}_0`$, since in products $`\stackrel{~}{G}_0\stackrel{~}{\mathrm{\Sigma }}\stackrel{~}{G}_0\mathrm{}\stackrel{~}{\mathrm{\Sigma }}\stackrel{~}{G}_0`$ only one simple pole at $`k=(0,\stackrel{~}{𝐤}_F)`$ survives, all other poles being cancelled by the corresponding zeros of the self-energy $`\stackrel{~}{\mathrm{\Sigma }}`$.

### B Superconducting state

To treat superconducting states we also add counterterms containing Cooper pair creation and annihilation operators, in addition to a shift of $`\xi _𝐤`$. We consider only spin singlet pairing, but triplet pairing could be dealt with analogously. We thus expand around a BCS mean-field Hamiltonian

$$\stackrel{~}{H}_0=\underset{𝐤,\sigma }{}\stackrel{~}{\xi }_𝐤n_{𝐤\sigma }+\underset{𝐤}{}\left[\mathrm{\Delta }_𝐤a_𝐤^{}a_𝐤^{}+\mathrm{\Delta }_𝐤^{}a_𝐤a_𝐤\right],$$
(7)
where $`\mathrm{\Delta }_𝐤`$ is the gap-function, which has to be determined self-consistently. In terms of Nambu operators

$$\mathrm{\Psi }_𝐤=\left(\begin{array}{c}a_𝐤\hfill \\ a_𝐤^{}\hfill \end{array}\right)\text{and}\mathrm{\Psi }_𝐤^{}=(a_𝐤^{},a_𝐤)$$
(8)
one can rewrite $`\stackrel{~}{H}_0`$ in more compact form as

$$\stackrel{~}{H}_0=\underset{𝐤}{}\stackrel{~}{\xi }_𝐤\mathrm{\Psi }_𝐤^{}\sigma _3\mathrm{\Psi }_𝐤\underset{𝐤}{}\mathrm{\Psi }_𝐤^{}(\mathrm{\Delta }_𝐤^{}\sigma _1\mathrm{\Delta }_𝐤^{\prime \prime }\sigma _2)\mathrm{\Psi }_𝐤,$$
(9)
where $`\sigma _1,\sigma _2,\sigma _3`$ are the Pauli matrices, and $`\mathrm{\Delta }_𝐤^{}`$ ($`\mathrm{\Delta }_𝐤^{\prime \prime }`$) is the real (imaginary) part of $`\mathrm{\Delta }_𝐤`$. The two expressions (7) and (9) for $`\stackrel{~}{H}_0`$ differ by the constant (c-number) $`_𝐤\stackrel{~}{\xi }_𝐤`$, which must be taken into account only when absolute energies are computed. The bare Nambu matrix propagator $`\stackrel{~}{𝐆}_0=\mathrm{\Psi }\mathrm{\Psi }^{}_{\stackrel{~}{0}}`$ following from $`\stackrel{~}{H}_0`$ is given by

$$\stackrel{~}{𝐆}_0^1(k)=\left(\begin{array}{cc}i\omega \stackrel{~}{\xi }_𝐤& \mathrm{\Delta }_𝐤\\ \mathrm{\Delta }_𝐤^{}& i\omega +\stackrel{~}{\xi }_𝐤\end{array}\right).$$
(10)
Extending the self-consistency condition for the normal state, we now require that the matrix self-energy $`\stackrel{\mathbf{~}}{𝚺}=\stackrel{~}{𝐆}_0^1𝐆^1`$ vanishes on the Fermi surface (defined by $`\stackrel{~}{\xi }_𝐤=0`$), that is

$$\stackrel{\mathbf{~}}{𝚺}(0,𝐤)=0\text{for}𝐤\stackrel{~}{}.$$
(11)
Thus, for $`\omega =0`$ and $`𝐤`$ on the Fermi surface, neither the diagonal nor the off-diagonal elements of $`\stackrel{~}{𝐆}_0^1(k)`$ are shifted by the interaction term $`\stackrel{~}{H}_I`$. The Feynman diagrams in Fig. 1 apply also to the superconducting case, if lines are interpreted as Nambu matrix propagators.

The above renormalized perturbation theory is reminiscent of the perturbation theory for symmetry broken phases formulated by Georges and Yedidia, where an order-parameter-dependent free energy function is constructed by adding Onsager reaction terms to the mean field contributions, and the actual order parameter is determined by minimizing this free energy.

## III Application to the Hubbard model

In this section we derive explicit expressions for the self-energy and the counterterms for the ground state $`(T=0)`$ of the Hubbard model, up to second order in the renormalized perturbation expansion.

The one-band Hubbard model

$$H=\underset{𝐢,𝐣}{}\underset{\sigma }{}t_{\mathrm{𝐢𝐣}}c_{𝐢\sigma }^{}c_{𝐣\sigma }+U\underset{𝐣}{}n_𝐣n_𝐣\mu N$$
(12)
describes lattice electrons with a hopping amplitude $`t_{\mathrm{𝐢𝐣}}`$ and a local interaction $`U`$. Here $`c_{𝐢\sigma }^{}`$ and $`c_{𝐢\sigma }`$ are creation and annihilation operators for electrons with spin projection $`\sigma `$ on a lattice site $`𝐢`$, and $`n_{𝐣\sigma }=c_{𝐢\sigma }^{}c_{𝐢\sigma }`$. Note that we have included the term $`\mu N`$ with the total particle number operator $`N`$ in our definition of $`H`$. The non-interacting part of $`H`$ can be written in momentum space as $`H_0=_𝐤\xi _𝐤n_{𝐤\sigma }`$ where $`\xi _𝐤=ϵ_𝐤\mu `$ and $`ϵ_𝐤`$ is the Fourier transform of $`t_{\mathrm{𝐢𝐣}}`$.

Our numerical results will be given for the Hubbard model on a square lattice with a hopping amplitude $`t`$ between nearest neighbors and a much smaller amplitude $`t^{}`$ between next-nearest neighbors. The corresponding dispersion relation is

$$ϵ_𝐤=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)4t^{}\mathrm{cos}k_x\mathrm{cos}k_y.$$
(13)
We now derive expressions for the self-energy and the resulting self-consistency equations up to second order in $`U`$.

### A Normal state

#### 1 First order

To first order in $`U`$ the self-energy is obtained as

$$\stackrel{~}{\mathrm{\Sigma }}^{(1)}(k)=U_k^{}\stackrel{~}{G}_0(k^{})e^{i\omega ^{}0^+}\delta \xi _𝐤,$$
(14)
where $`_k`$ is a short-hand notation for the frequency and momentum integral, including the usual factors of $`(2\pi )^1`$ for each integration variable. The first term results from diagram (1a) in Fig. 1, the second from diagram (1b). Note that the tadpole diagram (1a) yields a k-independent contribution, since the Hubbard interaction is local. The self-consistency condition (6) for $`\stackrel{~}{\mathrm{\Sigma }}^{(1)}`$ yields, after carrying out the $`\omega ^{}`$-integration,

$$\delta \xi _𝐤=U\frac{d^dk^{}}{(2\pi )^d}\mathrm{\Theta }(\mu ϵ_𝐤^{}\delta \xi _𝐤^{}),$$
(15)
to be satisfied (at least) for $`𝐤\stackrel{~}{}`$. Since the right hand side of this condition is a constant, it is natural to define $`\delta \xi _𝐤`$ by this constant for all $`𝐤`$. Using Luttinger’s theorem one can identify the above momentum integral with the particle density per spin, such that $`\delta \xi _𝐤=Un/2`$, where $`n`$ is the total density. The self-consistency condition thus yields the $`n(\mu )`$ relation of the interacting system. Since the counterterm can be chosen k-independently at first order, it may be interpreted as a shift of the chemical potential.

#### 2 Second order

The diagrams (2b) and (2c) from Fig. 1 obviously cancel each other to the extent that the first order diagrams (1a) and (1b) cancel. Writing $`\delta \xi _𝐤=\delta \xi ^{(1)}+\delta \xi _𝐤^{(2)}`$ with $`\delta \xi ^{(1)}`$ given by the constant on the right hand side of Eq. (15), such that $`\delta \xi _𝐤^{(2)}`$ is of order $`U^2`$ for all $`𝐤`$, the sum of contributions from (2b) and (2c) is of order $`U^3`$ and can thus be ignored at second order. Hence, only diagram (2a) contributes to the second order self-energy. Using the Feynman rules  one obtains

$$\stackrel{~}{\mathrm{\Sigma }}^{(2)}(k)=U^2_q\stackrel{~}{\mathrm{\Pi }}_0(q)\stackrel{~}{G}_0(kq),$$
(16)
where $`\stackrel{~}{\mathrm{\Pi }}_0(q)=_k^{}\stackrel{~}{G}_0(k^{})\stackrel{~}{G}_0(k^{}+q)`$. Adding first and second order terms, one arrives at the second order self-consistency condition

$$\delta \xi _𝐤=U\frac{d^dk^{}}{(2\pi )^d}\mathrm{\Theta }(\stackrel{~}{\xi }_𝐤^{})+\stackrel{~}{\mathrm{\Sigma }}^{(2)}(0,𝐤).$$
(17)
The counterterm $`\delta \xi _𝐤`$ has to be chosen such that the above equation is satisfied for all $`𝐤\stackrel{~}{}`$, that is for all $`𝐤`$ satisfying $`\stackrel{~}{\xi }_𝐤=0`$. Since $`\stackrel{~}{\mathrm{\Sigma }}^{(2)}(0,𝐤)`$ is momentum dependent, $`\delta \xi _𝐤`$ cannot be chosen constant any more. As a consequence, the Fermi surface of the interacting system will be deformed by interactions, even if the volume of the Fermi sea is kept fixed. Luttinger’s theorem can be used to determine the density from the volume of the Fermi sea as $`n=2\frac{d^dk}{(2\pi )^d}\mathrm{\Theta }(\stackrel{~}{\xi }_𝐤)`$.

### B Superconducting state

For the matrix elements of the Nambu propagator we use the standard notation

$$𝐆(k)=\left(\begin{array}{cc}G(k)& F(k)\\ F^{}(k)& G(k)\end{array}\right),$$
(18)
and the analogous expression for $`\stackrel{~}{𝐆}_0(k)`$. The matrix elements of the self-energy are denoted by

$$\stackrel{\mathbf{~}}{𝚺}(k)=\left(\begin{array}{cc}\stackrel{~}{\mathrm{\Sigma }}(k)& \stackrel{~}{S}(k)\\ \stackrel{~}{S}^{}(k)& \stackrel{~}{\mathrm{\Sigma }}(k)\end{array}\right).$$
(19)
#### 1 First order

In the presence of an off-diagonal counterterm $`\mathrm{\Delta }_𝐤`$, the diagonal part of $`\stackrel{\mathbf{~}}{𝚺}`$ is still given by Eq. (14) to first order, where $`\stackrel{~}{G}_0(k)`$ now depends on the gap function:

$$\stackrel{~}{G}_0(k)=\frac{i\omega +\stackrel{~}{\xi }_𝐤}{\omega ^2+\stackrel{~}{\xi }_𝐤^2+|\mathrm{\Delta }_𝐤|^2}.$$
(20)
The first order self-consistency relation (15) thus generalizes to

$$\delta \xi _𝐤=U\frac{d^dk^{}}{(2\pi )^d}\frac{1}{2}\left(1\stackrel{~}{\xi }_𝐤^{}/E_𝐤^{}\right),$$
(21)
with $`E_𝐤=\sqrt{\stackrel{~}{\xi }_𝐤^2+|\mathrm{\Delta }_𝐤^2|}`$. Note that the above integral is the BCS formula for the average particle density per spin.

The off-diagonal matrix element of $`\stackrel{\mathbf{~}}{𝚺}`$ is obtained from diagrams (1a) and (1b) in Fig. 1 as

$$\stackrel{~}{S}^{(1)}(k)=U_k^{}\stackrel{~}{F}_0(k^{})+\mathrm{\Delta }_𝐤$$
(22)
to first order in $`U`$, with

$$\stackrel{~}{F}_0(k)=\frac{\mathrm{\Delta }_𝐤}{\omega ^2+\stackrel{~}{\xi }_𝐤^2+|\mathrm{\Delta }_𝐤|^2}.$$
(23)
The off-diagonal part of the self-consistency condition (11) follows as

$$\mathrm{\Delta }_𝐤=U\frac{d^dk^{}}{(2\pi )^d}\frac{\mathrm{\Delta }_𝐤^{}}{2E_𝐤^{}}.$$
(24)
Extended as a condition for all $`𝐤`$ (and not just on $`\stackrel{~}{}`$) this is nothing but the gap equation for the Hubbard model as obtained by standard BCS theory. The self-consistency relation requires that $`\mathrm{\Delta }_𝐤`$ be constant on the Fermi surface, such that one naturally chooses a constant $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }`$ as an ansatz for all $`𝐤`$. A non-trivial solution $`\mathrm{\Delta }0`$ of this gap equation can obviously be obtained only for the attractive Hubbard model ($`U<0`$).

#### 2 Second order

The diagrams (2b) and (2c) cancel each other for the same reason as in the normal state. The contribution from diagram (2a) to the diagonal part of the self-energy is still given by formula (16), with $`\stackrel{~}{G}_0`$ from Eq. (20) and

$$\stackrel{~}{\mathrm{\Pi }}_0(q)=_k^{}\left[\stackrel{~}{G}_0(k^{})\stackrel{~}{G}_0(k^{}+q)+\stackrel{~}{F}_0(k^{})\stackrel{~}{F}_0^{}(k^{}+q)\right].$$
(25)
The second order contribution to the off-diagonal matrix element of $`\stackrel{\mathbf{~}}{𝚺}`$ is

$$\stackrel{~}{S}^{(2)}(k)=U^2_q\stackrel{~}{\mathrm{\Pi }}_0(q)\stackrel{~}{F}_0(kq).$$
(26)
The self-consistency relations read:

$`\delta \xi _𝐤`$ $`=`$ $`U{\displaystyle \frac{d^dk^{}}{(2\pi )^d}\frac{1}{2}\left(1\stackrel{~}{\xi }_𝐤^{}/E_𝐤^{}\right)}+\stackrel{~}{\mathrm{\Sigma }}^{(2)}(0,𝐤),`$ (27)
$`\mathrm{\Delta }_𝐤`$ $`=`$ $`U{\displaystyle \frac{d^dk^{}}{(2\pi )^d}\frac{\mathrm{\Delta }_𝐤^{}}{2E_𝐤^{}}}\stackrel{~}{S}^{(2)}(0,𝐤).`$ (28)
In Appendix A we present more explicit expressions for $`\stackrel{~}{\mathrm{\Sigma }}^{(2)}(0,𝐤)`$ and $`\stackrel{~}{S}^{(2)}(0,𝐤)`$, obtained by carrying out the frequency integrals.

### C Numerical solution

The self-consistency conditions are non-linear equations for the counterterms $`\delta \xi _𝐤`$ and, in the superconducting state, $`\mathrm{\Delta }_𝐤`$. The Fermi surface of the interacting system, $`\stackrel{~}{}`$, on which the self-consistency conditions must be satisfied, is not known a priori. The equations involve one momentum integral at first order, and two momentum integrals at second order. Such a non-linear system can only be solved iteratively. In this subsection we describe some details of our algorithm.

Since the counterterms are determined by the self-consistency conditions only on the Fermi surface, their momentum dependence away from $`\stackrel{~}{}`$ can be parametrized in many ways. We have chosen $`\delta \xi _𝐤`$ and $`\mathrm{\Delta }_𝐤`$ as constant along the straight lines connecting the square shaped line defined by the condition $`|k_x|+|k_y|=\pi `$ with the points $`(0,0)`$ and $`(\pi ,\pi )`$ of the Brillouin zone, respectively (see Fig. 2). For a numerical solution the remaining tangential momentum dependence is discretized by up to 256 points.

The iteration procedure starts with a tentative choice of counterterms. To be able to reach a symmetry broken solution one usually has to offer at least a small symmetry breaking counterterm in the beginning. In each iteration step new counterterms are determined via Eq. (17) in the normal state, and by Eqs. (27) and (28) for the superconducting state. The right hand side of these equations is evaluated using the counterterms obtained in the previous step, and $`𝐤`$ is chosen on the Fermi surface defined by the previous $`\delta \xi _𝐤`$. The momentum integrals are carried out using a Monte-Carlo routine. The iteration is continued until convergence is achieved, that is until the counterterms remain invariant within numerical accuracy from step to step. In all cases studied different choices of initial counterterms lead to the same unique solution. The symmetry breaking terms are much larger than the stochastic noise from the Monte-Carlo routine in all results shown.

The density is kept fixed by adjusting the chemical potential during the iteration procedure. To avoid a higher numerical effort we have computed the density from the Fermi surface volume in the normal state (justified by Luttinger’s theorem), and from the BCS formula for the density in superconducting solutions. The latter reduces to the Fermi surface volume in the normal state limit, such that the potential error of this approximation is very small as long as the gap is small.

## IV Results

We now discuss the most interesting results obtained within the renormalized perturbation theory described above, focussing mainly on the repulsive Hubbard model $`(U>0`$), for which we have found superconducting solutions with d-wave symmetry, as well as symmetry-breaking Fermi surface deformations.

### A Repulsive Hubbard model

The following results for the repulsive Hubbard model have been computed for the parameters $`t^{}=0.15t`$ and $`U=3t`$. The interaction is thus in the weak to intermediate coupling regime. For too small $`U`$-values it becomes very hard to resolve the small superconducting gap in the numerical solution.

We have solved the self-consistency equations for various densities ranging from $`n=0.88`$ to $`n=0.90`$, for which the Fermi surfaces are quite close to the saddle points of the bare dispersion relation $`ϵ_𝐤`$, located at $`(\pi ,0)`$ and $`(0,\pi )`$. In all cases the normal state is unstable towards superconductivity. The gap function in the superconducting state obtained from the self-consistency equations has $`d_{x^2y^2}`$-wave shape, with slight deviations from perfect d-wave symmetry in cases where the Fermi surface breaks the symmetry of the square lattice. This is in agreement with widespread expectations for the Hubbard model, and in particular with recent renormalization group arguments and calculations. In Fig. 3 we show the gap functions obtained at the densities $`n=0.88`$ and $`n=0.9`$, respectively. We note that the size of the gap is roughly one order of magnitude smaller than the critical cutoff scale $`\mathrm{\Lambda }_c`$ at which Cooper pair susceptibilities diverge in 1-loop renormalization group calculations for comparable model parameters. There are various possible reasons for this quantitative discrepancy. First, and probably most importantly, the enhancement of effective interactions due to fluctuations, especially antiferromagnetic spin fluctuations, is captured much better by a renormalization group calculation. Second, the approximate Fermi surface projection of vertices driving the renormalization group flow can lead to an overestimation of effective interactions and hence of critical energy scales. Furthermore, a renormalization group calculation within the symmetry broken phase could yield a gap that is somewhat smaller than $`\mathrm{\Lambda }_c`$.

While superconductivity is the only possible instability of the normal Fermi liquid state in the weak coupling limit (except for the case of perfect nesting at half-filling), at higher $`U`$ one should also consider the possibility of other, in particular magnetic, instabilities. This could be done within renormalized perturbation theory by allowing for counterterms introducing magnetic or charge order.

The Fermi surface is always deformed by interactions. The shifts generated by the momentum dependence of the counterterm $`\delta \xi _𝐤`$ are not very large. They are more pronounced near the saddle points of $`ϵ_𝐤`$, where small energy shifts lead to relatively large shifts in k-space. However, the results presented in Fig. 4 show that the Fermi surface of the interacting system can nevertheless differ strikingly from the bare one. For the densities $`n=0.880.889`$ the Fermi surface of the interacting system obviously breaks the point group symmetry of the square lattice. For $`n=0.88`$ and $`n=0.888`$ even the topology of the Fermi surface is changed by interactions. The deformed surface has open topology in these cases, instead of being closed around the points $`(0,0)`$ or $`(\pi ,\pi )`$ in the Brillouin zone. Note that the symmetry-broken Fermi surfaces shown here correspond to stable solutions of the self-consistency equations for the counterterms, while symmetric solutions are unstable.

More details about the Fermi surface shifts can be extracted from a plot of the second order counterterms, shown in Fig. 5. The actual shifts are determined by these terms plus a constant due to the first order counterterm and a shift of the chemical potential. At fixed density the interaction shifts the Fermi surface outwards at points where $`\stackrel{~}{\mathrm{\Sigma }}(0,\stackrel{~}{𝐤}_F)`$ has an absolute minimum, and inwards at points corresponding to absolute maxima. Interactions thus reduce the curvature of the Fermi surface near the diagonals in the Brillouin zone. Fig. 5 reveals that the Fermi surface deformation is slightly asymmetric also for $`n=0.9`$, but the symmetry breaking is too small to be seen in Fig. 4.

If the Fermi surface breaks the square lattice symmetry, the gap function $`\mathrm{\Delta }_𝐤`$ cannot have pure d-wave symmetry any more. See, for example, the gap function at density $`n=0.88`$ in Fig. 3. The deviation from perfect d-wave form is however quite small, since the symmetry breaking Fermi surface deformation is small.

Interaction-induced Fermi surface deformations which break the symmetry of the square lattice have already been discussed earlier in the literature. Yamase and Kohno  have obtained symmetry-broken Fermi surfaces within a slave boson mean-field theory for the t-J model. The effective interactions obtained from 1-loop renormalization group flows for the Hubbard model also favor symmetry-breaking Pomeranchuk instabilities of the Fermi surface, if the latter is close to the van Hove points. A systematic stability analysis of the Hubbard model using Wegner’s Hamiltonian flow equation method confirmed that symmetry breaking Fermi surface deformations are among the strongest instabilities. It remained an open question, however, whether such Fermi surface instabilities would be cut off by the superconducting gap. We have observed within our renormalized perturbation theory that symmetry breaking Fermi surface deformations occur indeed more easily, if the system is forced to stay in a normal state, by setting $`\mathrm{\Delta }_𝐤=0`$. Whether a symmetry broken Fermi surface and superconductivity coexist can be seen only by performing a calculation within the symmetry-broken state. This has not yet been done using the renormalization group or flow equation methods.

From a pure symmetry-group point of view the symmetry breaking generated by the Pomeranchuk instability is equivalent to that in ”nematic” electron liquids, first discussed by Kivelson et al.. These authors considered doped Mott insulators, that is strongly interacting systems. A general theory of orientational symmetry-breaking in fully isotropic (not lattice) two- and three-dimensional Fermi liquids has been reported by Oganesyan et al. Superconducting nematic states, in which discrete orientational symmetry breaking develops in addition to d-wave superconductivity, have been considered recently by Vojta et al. Motivated by experimental properties of single-particle excitations in cuprate superconductors they performed a general classification and field-theoretic analysis of various phases with an additional order parameter on top of $`d_{x^2y^2}`$-pairing.

### B Attractive Hubbard model

For the attractive Hubbard model ($`U<0`$) the renormalized perturbation expansion yields s-wave superconductivity already at first order, which is equivalent to BCS mean-field theory. At this level the gap function is constant in k-space. Extending the calculation to second order, a weak momentum dependence of $`\mathrm{\Delta }_𝐤`$ is generated, as seen in Fig. 6 for the parameters $`U=2t`$, $`t^{}=0.15t`$ and $`n=0.9`$. More importantly, the overall size of the gap is strongly reduced by fluctuations included in the second order terms. The average gap in Fig. 6 is only one third of the corresponding mean-field gap. It has been pointed out previously that fluctuations not contained in mean-field theory reduce the size of magnetic and other order parameters even in the weak coupling limit.

## V Conclusion

In summary, we have formulated a renormalized perturbation theory for interacting Fermi systems, which treats Fermi surface deformations and superconductivity via additive counterterms. This method is very convenient for studying the role of fluctuations for spontaneous symmetry breaking in a controlled weak-coupling expansion. A concrete application of the expansion carried out to second order yields several non-trivial results for the two-dimensional Hubbard model. In particular, for the repulsive model we have obtained the gap function of the expected d-wave superconducting state and, for Fermi levels close to the van Hove energy, an interacting Fermi surface with broken lattice symmetry, and in some cases even open topology. The symmetry-breaking pattern of the states with symmetry-broken Fermi surfaces is equivalent to that of ”nematic” electron liquids discussed already earlier from a different point of view.

The present work can be extended in several interesting directions. After fixing the counterterms one can compute the full momentum and energy dependence of the self-energy, and hence the spectral function for single-particle excitations. At second order the combined effects of symmetry breaking and quasi-particle decay are captured. Allowing for other symmetry-breaking counterterms, for example spin density waves, one can study the competition of magnetic, charge, and superconducting instabilities, as well as their possible coexistence. Finally, the formalism can be extended to finite temperature. In that case the singularities of the bare propagator are cut off by the smallest Matsubara frequency, but Fermi surface shifts and symmetry breaking can still be conveniently taken into account by counterterms.

Acknowledgements:
We are grateful to A. Georges, M. Keller, D. Rohe and M. Salmhofer for valuable discussions, and to D. Rohe also for a critical reading of the manuscript.

## A Frequency integrals

The Matsubara frequency integrals in the second order self-energy contributions can be carried out analytically by using the residue theorem. We only present the results for the superconducting case; the normal state results can be recovered by setting $`\mathrm{\Delta }_𝐤=0`$ in the following expressions.

The frequency integrals relevant for the evaluation of $`\stackrel{~}{\mathrm{\Pi }}_0`$ defined by Eq. (25) are

$`{\displaystyle \frac{dk_0}{2\pi }\stackrel{~}{G}_0(k)\stackrel{~}{G}_0(k+q)}`$ $`=`$ $`{\displaystyle \frac{E_𝐤+E_{𝐤+𝐪}}{2E_𝐤E_{𝐤+𝐪}}}{\displaystyle \frac{\stackrel{~}{\xi }_𝐤\stackrel{~}{\xi }_{𝐤+𝐪}E_𝐤E_{𝐤+𝐪}}{q_0^2+[E_𝐤+E_{𝐤+𝐪}]^2}}`$ (A1)
$`+`$ $`{\displaystyle \frac{iq_0}{2E_𝐤E_{𝐤+𝐪}}}{\displaystyle \frac{\stackrel{~}{\xi }_𝐤E_{𝐤+𝐪}E_𝐤\stackrel{~}{\xi }_{𝐤+𝐪}}{q_0^2+[E_𝐤+E_{𝐤+𝐪}]^2}}`$ (A2)
and

$$\frac{dk_0}{2\pi }\stackrel{~}{F}_0(k)\stackrel{~}{F}_0^{}(k+q)=\frac{E_𝐤+E_{𝐤+𝐪}}{2E_𝐤E_{𝐤+𝐪}}\frac{\mathrm{\Delta }_𝐤\mathrm{\Delta }_{𝐤+𝐪}^{}}{q_0^2+[E_𝐤+E_{𝐤+𝐪}]^2}.$$
(A3)
The imaginary part of $`\stackrel{~}{\mathrm{\Pi }}_0`$ does not contribute to $`\stackrel{\mathbf{~}}{𝚺}(0,𝐤)`$. Carrying out the $`q_0`$-integral in Eqs. (16) and (26) yields:

$`\stackrel{~}{\mathrm{\Sigma }}^{(2)}(0,𝐤)`$ $`=`$ $`U^2{\displaystyle _𝐪}{\displaystyle _𝐤^{}}\stackrel{~}{\xi }_{𝐤𝐪}C(𝐤,𝐤^{},𝐪)`$ (A4)
$`\stackrel{~}{S}^{(2)}(0,𝐤)`$ $`=`$ $`U^2{\displaystyle _𝐪}{\displaystyle _𝐤^{}}\mathrm{\Delta }_{𝐤𝐪}C(𝐤,𝐤^{},𝐪),`$ (A5)
where

$$C(𝐤,𝐤^{},𝐪)=\frac{E_𝐤^{}E_{𝐤^{}+𝐪}\stackrel{~}{\xi }_𝐤^{}\stackrel{~}{\xi }_{𝐤^{}+𝐪}\mathrm{\Delta }_𝐤^{}\mathrm{\Delta }_{𝐤^{}+𝐪}^{}}{4E_{𝐤𝐪}E_𝐤^{}E_{𝐤^{}+𝐪}\left[E_{𝐤𝐪}+E_𝐤^{}+E_{𝐤^{}+𝐪}\right]}.$$
(A6)
FIGURES