text:
# Giant fluctuations of superconducting order parameter in ferromagnet/superconductor single electron transistors

## Abstract

Spin dependent transport in a ferromagnet/superconductor/ferromagnet single electron transistor is studied theoretically with spin accumulation, spin relaxation, gap suppression, and charging effects taken into account. A strong dependence of the gap on the magnetic state of the outer electrodes is found, which gives rise to a negative magneto-resistance of up to 100 %. We predict that fluctuations of the spin accumulation due to tunneling of quasi-particles can play such an important role as to cause the island to fluctuate between the superconducting and normal state.

Electronic devices, which exploit the spin as well as charge of the electron have attracted much interest over the last several years. This interest in “spintronics” is stimulated by expected industrial applications, such as magnetic tunnel junctions based sensors, magnetic random access memory, as well as basic physics, offering new ways to manipulate quantum states using these two fundamental properties of the electron.

In nanoscale magnetic devices, such as single electron transistors (SET’s), the presence of Coulomb blockade makes the charge transport sensitive to a very few electron states, which results in a number of pronounced effects in magneto-conductance (see review Martinek:jap:review,  and refs. therein). The sensitivity to spin imbalance in the SET should be greatly enhanced if the island is made superconducting, since the dependence of the gap on the number of accumulated spins is highly non-linear. Takahashi *et al.* Takahashi  analysed a double tunnel barrier structure with ferromagnetic outer layers and superconducting central layer (F/S/F) without charging effects. They found that the anti-parallel alignment of the ferromagnetic films can produce a non-equilibrium spin distribution in the central layer, which suppresses the superconducting gap and thereby strongly influences the transport. Immamura *et al.* Imamura\_FSF\_parity  considered the parity effect in a ferromagnetic SET with a superconducting island (F/S/F SET) in the limiting case of charging energy larger than the superconducting gap, low voltage, and zero temperature. They found that the spin accumulation is zero in the limit considered, i.e., the superconducting gap of the island is unaffected. Magnetoresistance (MR) can exist due to the polarized spin of the single excess electron but is reduced in magnitude from the Julliere value by spin relaxation.

In this letter we develop a comprehensive model of the F/S/F SET by taking into account spin accumulation, spin relaxation, gap suppression, charging effects as well as conduction fluctuations in experimentally realistic structures, at finite temperatures. We find that spin accumulation in the F/S/F SET suppresses the gap, which strongly affects the transport properties and results in a greatly enhanced MR. Progressively smaller spin accumulation is required to suppress the gap as the size of the island is reduced, making it feasible to electrically detect just a few non-equilibrium spins. We find that fluctuations of the spin accumulation can play such an important role as to cause giant fluctuations of the gap between the full value at a given temperature and zero (between the superconducting and normal state *at fixed bias and gate voltages*), with corresponding large fluctuations in the conductance of the device. Even in the parallel magnetic state of the device, where the average spin accumulation is zero, fluctuations of spin accumulation can lead to a significant gap suppression. The transistor exhibits a practically full *gate-controlled* MR effect (MR$`1`$, corresponding to $`1000`$% in typical GMR notations). We compare the results of the model with the recent, rather divergent experimental findings on ferromagnetic SET’s Chen ; Johansson ; Johansson\_comment .

The model geometry is shown in Fig.1(a). Two ferromagnetic electrodes are connected to a superconducting island through tunnel junctions. The tunnel resistance is assumed to be much larger than the Quantum resistance, $`R_iR_Q=h/e^226\text{k}\mathrm{\Omega }`$. We assume collinear magnetization of the two electrodes, the magnetization of the left electrode can be either parallel (*P*) or anti-parallel (*AP*) to the magnetization of the right electrode. The resistance of the tunnel junctions is spin dependent: $`R_i^{}=2R_i/(1+P_i)`$, $`R_i^{}=2R_i/(1P_i)`$ (where $`i=L`$ (left) or $`R`$ (right) denotes the junction). Here the polarization is of the same sign ($`P_LP_R>0`$) for the parallel alignment. Letting $`P_LP_L`$ (or $`P_RP_R`$) for reversed magnetization of the left (or right) electrode, we obtain the resistances in the anti-parallel alignment.

The spin relaxation time can be comparable to or greater than the average time between tunnel events, so that a non-equilibrium spin distribution can accumulate on the island, resulting in a shift in the chemical potential. The shift in chemical potential is schematically shown in Fig.1(b). We consider the regime where the energy relaxation time is much shorter than the average time between two successive tunnel events (no gap suppression is expected in the opposite, elastic regime; for details see Tserkovnyak:FSF, ). The electrons are then instantaneously thermalized after they have tunneled and are described by the Fermi distribution function. This shift in chemical potential, $`\mathrm{\Delta }\mu ^\sigma `$, is related to the deviation in the number of spins from the equilibrium value, $`n^\sigma `$, for each spin-band, $`\sigma =,`$:

$$n^\sigma =\rho V_{\mathrm{}}^{\mathrm{}}𝑑EN(E)\left[f(E\mathrm{\Delta }\mu ^\sigma )f(E)\right],$$
(1)
where $`N(E)=\mathrm{\Theta }(E\mathrm{\Delta })\frac{|E|}{\sqrt{E^2\mathrm{\Delta }^2}}`$ is the BCS density of states, $`\rho `$ is the density of states at the Fermi energy (per spin) in the normal state, $`f(E)`$ is the Fermi distribution function, $`V`$ is the island volume. The shift in the chemical potential on the superconducting island influences the gap through the BCS gap equation Tinkham

$$\mathrm{ln}\left(\frac{\mathrm{\Delta }_0}{\mathrm{\Delta }}\right)=_0^{\mathrm{}\omega _D}\frac{dE}{\xi }\left[f(\xi +\mathrm{\Delta }\mu ^{})+f(\xi +\mathrm{\Delta }\mu ^{})\right],$$
(2)
where $`\xi =\sqrt{E^2+\mathrm{\Delta }^2}`$ and $`\omega _D`$ is the Debye frequency. For arbitrary spin distribution on the island, (arbitrary $`n^{}`$ and $`n^{}`$), solving equations (1) and (2) self-consistently yields the shift in the chemical potential, $`\mathrm{\Delta }\mu ^\sigma `$ , and the gap, $`\mathrm{\Delta }`$. In Fig.2 the chemical potential and the gap as a function of the spin imbalance, $`sn^{}n^{}`$, are shown for the case $`n^{}=n^{}`$. When the spin imbalance increases, there will be a critical value of $`s`$ where the gap will discontinuously drop to zero and the chemical potential will drop to a lower value. After the complete suppression of the gap the chemical potential will increase linearly with the slope of $`1/\rho V`$ as in the normal conducting case. Reducing the spin imbalance will restore the gap at the same critical value of $`s`$. Note that the gap and chemical potential do not only depend on $`n^{}`$ and $`n^{}`$ through the combination $`s=n^{}n^{}`$, the total spin accumulation. Different $`n^{}`$ and $`n^{}`$ having the same $`s`$ can produce different $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }\mu ^\sigma `$.

The transport properties are determined by the tunneling rates according to the “orthodox theory” SCT . The tunneling rates depend on the charge of the island through $`n=n^{}+n^{}`$, as well as $`n^{}`$ and $`n^{}`$ separately, through the dependence of $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }\mu ^\sigma `$ on $`n^{}`$ and $`n^{}`$. The tunneling rates are spin dependent also through the spin dependent resistances:

$`\mathrm{\Gamma }_i^{\text{on},\sigma }(n)`$ $`=`$ $`{\displaystyle \frac{1}{e^2R_i^\sigma }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑EN(E)f(E\delta E_i^{\text{on}}(n))\left[1f(E+\mathrm{\Delta }\mu ^\sigma )\right]`$
$`\mathrm{\Gamma }_i^{\text{off},\sigma }(n)`$ $`=`$ $`{\displaystyle \frac{1}{e^2R_i^\sigma }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑EN(E)f(E\mathrm{\Delta }\mu ^\sigma )\left[1f(E\delta E_i^{\text{off}}(n))\right],`$
where $`\delta E_i^{\text{on/off}}(n)`$ is the change in the charging energy for tunneling one electron on or off the island through the left or right junction:

$$\delta E_i^{\text{on/off}}(n)=\frac{e^2}{C_\mathrm{\Sigma }}\left[\pm \left(nn_g+\frac{V_i}{e}\frac{C_LC_R}{C_i}\right)+\frac{1}{2}\right]$$

Here $`C_\mathrm{\Sigma }=C_L+C_R+C_g`$ is the total capacitance of the island and $`n_g=C_gV_g/e`$ is the normalized gate charge. Tunneling that increases the the energy have exponentially small rates (zero at $`T=0`$), giving the characteristic Coulomb blockade below a threshold voltage.

In what follows we choose the parameter space to represent the typical nano-lithographically fabricated SET structures. We therefore consider SET’s with island dimensions in the $`10^510^6`$ nm3 range, such that up to a hundred states will participate in the transport. The relatively large number of states involved makes the Monte Carlo method Likharev  more suitable than the master equation approach. The intertwined dependence of the gap, the chemical potential and the spin accumulation can be included in a straight forward way with the Monte Carlo method. Moreover, the fundamental stochastic nature of tunneling and the resultant fluctuations in the spin transport through the structure are well modeled by the Monte Carlo procedure.

The state of the system is determined by $`n^{}`$ and $`n^{}`$, from which the gap and shift in chemical potential are calculated using Eqs. (1) and (2). After obtaining $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }\mu ^\sigma `$, the eight possible tunneling rates are computed, which increase or decrease $`n^{}`$ or $`n^{}`$ by one. The interval zero to one is then divided into eight segments with widths proportional to the tunneling rates and a random number in the interval  is generated to choose one of the eight tunnel events. The probability to preserve the present configuration decays as $`P=\mathrm{exp}[(tt_0)\mathrm{\Gamma }_\mathrm{\Sigma }]`$, where $`t_0`$ is the time of the previous tunnel event and $`\mathrm{\Gamma }_\mathrm{\Sigma }`$ is the sum of the eight tunneling rates. A random number in the interval  gives the time between successive tunnel events by inverting the equation for the probability. The spin flip process is then accounted for in a phenomenological way by letting $`ss\mathrm{exp}[(tt_0)/\tau _{\text{sf}}]`$, where $`\tau _{\text{sf}}`$ is the spin flip time. This procedure is repeated, keeping track of all quantities, until good statistics is obtained.

The dynamic effect of gap suppression is shown in Fig.3, where the spin accumulation and superconducting gap are plotted versus time at fixed bias and gate voltages. As soon as the spin imbalance reaches the critical value ($`18`$ for these device parameters) the gap drops to zero and the island goes normal. At the critical bias, fluctuations of only a few polarized spins on the island cause giant fluctuations of the order parameter from the full value at the given temperature to zero. This effect should be directly observable in experiments employing time resolved (\<1$`\mu `$s scale) conductance measurements. Typical current-voltage (I-V) measurements (\>1 ms time constant) will record the gap, which is a time average of the dynamic gap shown in Fig.3. Fig.4 shows this mean value of the superconducting gap, $`\mathrm{\Delta }`$, as well as its standard deviation, $`\delta \mathrm{\Delta }=[\mathrm{\Delta }^2\mathrm{\Delta }^2]^{1/2}`$, as a function of bias voltage. The fluctuations result in a smooth rather than discontinuous transition from the superconducting to the normal state as the spin accumulation is increased in the *AP* state by increasing the bias. The fluctuations peak at the critical bias in the *AP* state (0.42 mV in this case). Contrary to what would be intuitively expected, fluctuations of the spin accumulation on the island in the *P* state can be very large. The fluctuations increase with bias voltage and becomes as large as the critical spin accumulation giving a significant gap suppression. This unexpected effect is a consequence of the stochastic dynamics of *spin*-dependent transport through the device (charge fluctuation is suppressed by the charging energy) and the highly non-linear dependence of the gap on spin accumulation, since the time averaged spin accumulation in the parallel state is zero.

The suppression of the superconducting gap by spin imbalance in the *AP* state results in fluctuations of the conductance of the device, which are greater in magnitude as the conductance itself. Fig.5a shows the time averaged I-V characteristics of the ferromagnetic SET in both the *AP* and *P* states. At the critical bias, the gap suppression leads to a 3-fold increase in current as a function of the magnetic state of the device, which results in a correspondingly large negative magnetoresistance. An even stronger MR effect, MR\>0.8 is obtained at the critical bias by applying a gate voltage as shown in Fig.5b. The island size chosen for the illustration above is representative of typical experimental implementations today, but not the smallest fabricated to date. For smaller superconducting islands, progressively fewer spins are needed to produce the full effect of gap suppression, with MR reaching $`1`$. Increasing the spin relaxation time or reducing thermal excitations by going to lower temperatures enhances the effect.

To estimate the parameters of a realistic experimental realisation of a F/S/F SET, assume a SET fabricated with Cobalt leads having a spin polarization of 0.35, Aluminum as the superconducting island with $`\mathrm{\Delta }=200\mu `$eV, $`\rho =10^{28}`$ \[m3eV\]-1 and a volume of $`15\times 70\times 200`$ nm3, so that $`\rho V=2.110^6`$ \[eV\]-1. The spin relaxation time of quasi-particles in superconducting *Al* has not been measured, but for normal Aluminum at 4.2 K spin relaxation times between 10 ns Johnson  and recently 0.1 ns Jedema  have been reported. Assuming a spin relaxation time of $`\tau _{\text{sf}}`$ 1 ns, $`T=`$0.05 K, $`R=30`$ k$`\mathrm{\Omega }`$, and $`E_\text{c}=100`$ $`\mu `$eV this sample would have the transport characteristics shown in Figs.2-5. If the spin relaxation time was shorter, the island volume would have to be smaller to achieve the same effect. Increased $`\tau _{sf}`$ in the superconducting state has been predicted theoretically Yamashita , which would lead to stronger gap suppression and higher MR.

When we compare our theoretical results with the recent experiments on F/S/F SET’s Chen ; Johansson , we draw the conclusion that the reported MR is unlikely due to gap suppression by spin accumulation. The volume of the islands studied was rather large in both experiments ($`\rho V10^7`$ \[eV\]-1) and the temperature relatively high (T $``$ 0.25 K). The combined effect of large volume and high temperature would give insignificant gap suppression in the voltage interval where large MR was reported. There are several other points where our simulations do not agree with the interpretation of experiments reported. In both experiments the maximum MR was observed at or below the threshold voltage, $`eV_{\text{th}}2\mathrm{\Delta }+E_\text{c}`$. In our calculations we always find $`V_{\text{th}}`$ to be the same in the two states, and the maximum MR occurs above $`V_{\text{th}}`$. This arises form the fact that the spin accumulation is directly proportional to the current, so there is no mechanism for spin accumulation below the threshold voltage, where no current flows. Further, the extracted gap in Chen  did not go to zero monotonically with increasing voltage but had a maximum at an intermediate voltage, contrary to what we predict. The temperature dependence of the MR measured in Johansson , with the peak of the MR occurring at the same bias voltage, is in contrast with our calculations where the MR peak shifts to higher bias with increasing temperature (not shown). More experiments are needed to resolve these controversies.

###### Acknowledgements.

J. Johansson gratefully acknowledge support from the Swedish SSF Spintronics program.


format: text

text:
# Nucleon Decay with Domain-Wall Fermions

## 1 INTRODUCTION

There are two methods to calculate nucleon decay matrix elements using lattice QCD. One is the direct calculation of the matrix elements dealing with three- and various two-point functions . The other is the indirect method where one relies on chiral perturbation theory whose low energy constants are calculated on the lattice . The former method is preferable but requires huge computation, typically 10 times larger than the latter to achieve similar statistical accuracy.

In a previous work , we made use of good chiral symmetry of the domain-wall fermion with the DBW2 gauge action . The proton decay matrix element of the pion final state was calculated with both methods. The results are consistent with each other within the relatively large error. There we assumed the mixing of operators arising from the explicit chiral symmetry breaking is negligible. Here we reconsider this issue by calculating the mixing with a non-perturbative renormalization (NPR). Also the matrix elements of the proton decay in case of the non-degenerate quark mass in the pseudoscalar state are calculated allowing further discussion on the effectiveness of the chiral perturbation theory.

All of our calculations are done on the quenched DBW2 gauge configurations with size $`N_\sigma ^3\times N_\tau =16^3\times 32`$. The domain-wall height is $`M=1.8`$, the fifth dimension size is $`L_s=16`$ for the NPR and $`12`$ for the matrix element.

## 2 OPERATOR MIXING

The nucleon decay is induced by the dimension-six operator which consists of one lepton and three quarks. The three quark operator, with which we calculate matrix element with initial nucleon and final pseudoscalar state, in general takes a form as

$$𝒪_{uds}^a=ϵ^{ijk}(u^{iT}C\mathrm{\Gamma }d^j)\mathrm{\Gamma }^{}s^k.$$
(1)
The superscript $`a`$ specifies a combination of the gamma matrix $`(\mathrm{\Gamma }\mathrm{\Gamma }^{})`$ listed in Table 1. Operators are classified by the ordinary parity ($`𝒫`$) and by the parity of switching ($`𝒮`$) $`u`$ and $`d`$ quarks. The $`𝒮`$ is a symmetry of the (lattice) QCD Lagrangian when quark masses are degenerate. The $`u`$, $`d`$, $`s`$ do not necessarily mean the real flavor, but simply stand for the different flavors. The other arrangements of flavors, such like $`𝒪_{usd}^a`$ and $`𝒪_{dus}^a`$ are always expressed in terms of $`𝒪_{uds}^a`$ by Fierz rearrangement and/or by $`𝒮`$. Operators with two flavors like $`𝒪_{udu}`$ follows the same argument.

Lattice operators are renormalized as

$$𝒪_{ren}^a=Z^{ab}𝒪_{latt}^b.$$
(2)
For example, for operators in $`𝒫^+`$, $`𝒮^{}`$ sector, which will not mix the operators in the other sectors, $`a`$ and $`b`$ in eq. 2 may be $`SS`$, $`PP`$, $`AA`$. The same $`3\times 3`$ matrix $`Z`$ applies to the $`𝒫^{}`$, $`𝒮^{}`$ sector: $`SP`$, $`PS`$, $`AV`$, too.

For the NPR with RI-MOM scheme  we calculate the Green’s function of the operator with the three quark external states in the Landau gauge,

$$G^a(x_0,x_1,x_2,x_3)=𝒪_{uds}^a(x_0)\overline{u}(x_1)\overline{d}(x_2)\overline{s}(x_3).$$
(3)
We set $`x_0=0`$. Fourier transformation is performed on the other legs with the same momentum $`p`$, which are then amputated to obtain the vertex function,

$$\mathrm{\Lambda }^a(p^2)=\text{F.T.}G^a(0,x_1,x_2,x_3)|_{Amp}.$$
(4)
Writing the tensor indices explicitly, the renormalization condition of the RI-MOM scheme leads,

$$P_{ijk\beta \alpha \delta \gamma }^aZ_q^{3/2}Z^{bc}\mathrm{\Lambda }_{ijk\alpha \beta \gamma \delta }^c=\delta ^{ab},$$
(5)
where $`Z_q`$ is the quark wave function renormalization $`i,j,k`$ are color indices, $`\alpha `$, $`\beta `$ and $`\gamma `$, $`\delta `$ are Dirac indices associated with $`\mathrm{\Gamma }`$, $`\mathrm{\Gamma }^{}`$ respectively. The renormalization condition should be applied where the pion pole effect  is absent and the lattice artifact is small : $`\mathrm{\Lambda }_{QCD}|p|1/a`$. The projection matrix $`P^a`$, for example onto $`SS`$, is $`P^{SS}=ϵ^{ijk}(C^1)^{\beta \alpha }\delta ^{\delta \gamma }/(644)`$. We define the matrix $`M`$ as,

$$M^{ab}=P_{ijk\beta \alpha \delta \gamma }^a\mathrm{\Lambda }_{ijk\alpha \beta \gamma \delta }^b,$$
(6)
which will give $`Z_q^{3/2}(Z^1)^{ab}`$.

Fig. 1 shows the resulting $`M^{ab}`$ at a quark mass $`m_f=0.025`$ for the domain-wall fermion. 51 configurations have been analyzed .

The renormalization factor of $`SS,SS`$ and $`PP,PP`$ are identical within the error except for the vicinity of zero momentum, the others are zero within the two standard deviation level. Absence of the off-diagonal element and coincidence of the $`SS,SS`$ and $`PP,PP`$ element are expected in the case of exact chiral symmetry. We have examined four different masses $`m_f=0.025,0.04,0.055,0.07`$, where mass dependences are hardly seen in all range of momentum except near zero.

As a result, the nucleon decay operator $`𝒪_{R/L;L}^{B/}=ϵ^{ijk}(u^{iT}CP_{R/L}d^j)P_Ls^k`$ \[$`S\pm PSP`$\] is renormalized multiplicatively without any mixing at practical level for this parameter of the domain-wall fermion. In the next section we show the results of the bare matrix elements with $`L_s=12`$. The residual mass $`m_{res}=1.6`$ MeV at $`L_s=12`$ is still comparable to $`0.7`$ MeV at $`L_s=16`$ . Hence the effect of $`L_s`$ on the renormalization factor is expected to be absent in this case, too, which needs to be checked by further performing simulation with $`L_s=12`$.

## 3 MATRIX ELEMENTS

We follow the ref.  to calculate the nucleon decay matrix elements by the direct method. Our results for the decay to a meson with degenerate quark mass are already shown in , where we have used 100 gauge configurations. We extend our calculation to the non-degenerate quark mass case to treat the decay to kaon. Our masses are $`m_{f1},m_{f2}=0.02,0,04,0.06,0.08`$ with $`m_{f1}m_{f2}`$, where $`m_{f1}`$ is for the $`u`$, $`d`$ quarks, $`m_{f2}`$ for the $`u`$, $`d`$ or $`s`$ quark. The other parameters are the same as in . Result of the relevant form factor $`W_0`$ of the matrix element for various decay modes are shown in Fig. 2. We are assuming the SU$`(2)`$ symmetry for the $`u`$ and $`d`$ quarks. There are other possible matrix elements, but they can be calculated with the matrix elements in the figure when the SU$`(2)`$ symmetry is intact.

The indirect method uses the chiral perturbation theory, which relates the low energy parameters in the chiral Lagrangian to the matrix elements. JLQCD extended the calculation by Claudson et. al  to all possible decay modes . The low energy parameters are mass ratio of nucleon and baryon with strangeness $`S=1`$, $`m_N/m_B`$, pion decay constant $`f`$, $`D`$ and $`F`$ parameterizing the semileptonic decay of baryons, nucleon decay parameter $`\alpha `$ and $`\beta `$. These parameters are in principle calculable on the lattice. As we are using the chiral fermion, observables calculated on the lattice should obey the chiral perturbation whose parameters are calculated on the lattice. However we use the experimental values except for $`\alpha `$ and $`\beta `$. The reasons are that we know the results of $`f_\pi `$  and $`g_A`$  $`(=D+F)`$ are anyway consistent with the experiment, and that we have no estimate of $`m_B`$ and individual $`D`$ and $`F`$. Using the $`\alpha `$ and $`\beta `$ in our previous study  and $`D=0.8,F=0.47`$ , $`f=0.131`$ GeV, $`m_N/m_B=0.817`$, the matrix elements are calculated and shown in Fig. 2. The lattice cutoff calculated by $`\rho `$ mass input, $`a^1=1.31`$ GeV is used to make dimensionless value.

The indirect method gives quite good estimate of the matrix element. Consistency of the two method within the error is in contrast to the result of JLQCD  with Wilson fermion. However, this could be caused by larger statistical error in our calculation. We need more statistics to be able to make more reliable statements.

## 4 SUMMARY

We have investigated the mixing of three-quark nucleon decay operators with domain-wall fermions by NPR. The mixing is practically negligible at the parameter set similar to that where the matrix element are calculated. We have extended our nucleon decay matrix element to those for the final kaon state matrix elements. The direct and indirect methods give consistent results within our statistical error. Performing NPR at the parameter of matrix element calculation together with the matching calculation will give us the matrix element in $`\overline{\text{MS}}`$.


text:
# 1 Introduction

## 1 Introduction

Our present understanding of the Universe is based upon the successful hot Big Bang theory, which explains its evolution from the first fraction of a second to our present age, around 13 billion years later. This theory rests upon four strong pillars, a theoretical framework based on general relativity, as put forward by Albert Einstein  and Alexander A. Friedmann  in the 1920s, and three robust observational facts: First, the expansion of the Universe, discovered by Edwin P. Hubble  in the 1930s, as a recession of galaxies at a speed proportional to their distance from us. Second, the relative abundance of light elements, explained by George Gamow  in the 1940s, mainly that of helium, deuterium and lithium, which were cooked from the nuclear reactions that took place at around a second to a few minutes after the Big Bang, when the Universe was a few times hotter than the core of the Sun. Third, the cosmic microwave background (CMB), the afterglow of the Big Bang, discovered in 1965 by Arno A. Penzias and Robert W. Wilson  as a very isotropic black-body radiation at a temperature of about 3 degrees Kelvin, emitted when the Universe was cold enough to form neutral atoms and photons decoupled from matter, approximately 400,000 years after the Big Bang. Today, these observations are confirmed to within a few per cent accuracy, and have helped establish the paradigm of inflationary hot Big Bang cosmology as the preferred model of the Universe .

## 2 The Universe as we see it

All of modern cosmology stems essentially from an application of the Copernican principle: we are not at the centre of the Universe. In fact, today we take Copernicus’ idea one step further and assert the “cosmological principle”: nobody is at the centre of the Universe. The cosmos, viewed from any point, looks the same as when viewed from any other point. This, like other symmetry principles more directly familiar to particle physicists, turns out to be an immensely powerful idea. In particular, it leads to the apparently inescapable conclusion that the Universe has a finite age. There was a beginning of time.

This cosmological principle can be expressed by applying one of the most “reasonable” symmetries to the Einstein equation, which initially connects the geometry of space-time with the matter content of the Universe. The simplest symmetries of this type are homogeneity  and isotropy. By homogeneity, we mean that the Universe is invariant under spatial translations, while isotropy implies invariance under rotations. In this sense the contents of the Universe can be modelled as a perfect fluid with some equation of state 111The equation of state is frequently given by a relation between density $`\rho `$ and pressure $`p`$ of the fluid.. While this is certainly a poor description of the contents of the Universe on small scales, such as the size of people or planets or even galaxies, it is an excellent approximation if we average over extremely large scales in the Universe, for which the matter is known observationally to be very smoothly distributed.

General relativity combined with homogeneity and isotropy leads to a startling conclusion: space-time is dynamic. The Universe is not static, but is bound to be either expanding or contracting. Indeed, in 1929, Hubble undertook a project to measure the distances to the spiral “nebulae”, as they had been known. Hubble’s method involved using Cepheid variables. Cepheid variables have the useful property that the period of their variation is correlated to their absolute brightness. Therefore, by measuring the apparent brightness and the period of a distant Cepheid, one can determine its absolute brightness and therefore its distance. Hubble applied this method to a number of nearby galaxies, and determined that almost all of them were receding from the Earth. Moreover, the more distant the galaxy was, the faster it was receding, according to a roughly linear relation: $`v=H_0d.`$ This is the famous Hubble law, and the constant $`H_0`$ is known as Hubble’s constant. The current best estimate, determined using the Hubble space telescope to resolve Cepheids in galaxies at unprecedented distances, is $`H_0=71\pm 6\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ 222The parsec is a spatial astronomical unit corresponding to one second of arc of the parallax measured from opposite sides of the earth’s orbit: $`1\mathrm{pc}=310^{18}\mathrm{cm}`$. . In any case, the Hubble law is exactly what one would expect from the so-called expanding Friedmann–Robertson–Walker (FRW) Universe .

The expansion of the Universe leads to a number of interesting effects. One is the cosmological redshift of photons. The usual way to see this is that, from the Hubble law, distant objects appear to be receding at a velocity $`v=H_0d`$, which means that photons emitted from the body are redshifted by the recession velocity of the source. There is another way to look at the same effect: because of the expansion of space, the wavelength of a photon increases with the scale factor: $`\lambda a(t),`$ so that a photon, propagating in the space as the universe expands, gets shifted to longer and longer wavelengths. The redshift $`z`$ of a photon is then given by the ratio of the scale factor today to the scale factor when the photon was emitted: $`1+z=\frac{a(t_0)}{a(t_{\mathrm{em}})}.`$ This redshifting due to expansion applies to particles other than photons as well. For some massive body moving relative to the expansion with some momentum $`p`$, the momentum also “redshifts”: $`p\frac{1}{a(t)}.`$ We then have the remarkable result that freely moving bodies in an expanding Universe eventually come to rest relative to the expanding coordinate system. Thus the expansion of the Universe creates a kind of dynamical friction for everything moving in it. If we take a bunch of particles with temperature $`T`$ in thermal equilibrium, the momenta of all these particles will decrease linearly with the expansion, and the system will cool. For a gas in thermal equilibrium, the temperature is in fact inversely proportional to the scale factor: $`T\frac{1}{a(t)}.`$ The current temperature of the Universe is $`2.725\pm 0.005\mathrm{K}`$ .

One of the things that cosmologists most want to measure accurately is the total density $`\rho `$ of the Universe. This is very often expressed in units of the critical density needed to make the geometry of the Universe flat. Observers have made attempts to measure the density of the Universe using a variety of methods, including measuring galactic rotation curves, the velocities of galaxies orbiting in clusters, X-ray measurements of galaxy clusters, the velocities and spatial distribution of galaxies on large scales, and gravitational lensing. These measurements have repeatedly pointed to the existence of a large amount of dark matter, different from the baryonic constituents of stars and planets. Moreover two groupes - the Supernova Cosmology Project  and the High-z Supernova Search Team  has recently arrived at a picture of the Universe with flat geometry: matter, including both baryons and the mysterious dark matter, makes up only about 30% of the energy density in the Universe. The remaining 70% is made of something that looks very much like a fluid with negative pressure, which slightly accelerate the current expansion of the Universe 333Deeper and more precies supernova measurements  confirm this conclusion. Future data on supernova from space based  and ground based  telescopes promise to yield substantial information on the history of cosmic acceleration, and to constrain possible models for the dark energy.. This dark energy can possibly be identified with the vacuum energy predicted by quantum field theory, except that the energy density is 120 orders of magnitude smaller than would be expected from a naive analysis. This density budget of the Universe, where the critical density is distributed between the following contributions: $`4.7\pm 0.6`$% luminous matter, $`24\pm 7`$% dark matter and $`71.3\pm 8`$% of dark energy, gets confirmed at very high confidence level by recent fascinating results of the WMAP satellite, which in fact ”took a picture” of the Universe when it was only about 400,000 years old.

## 3 Thermal history of the Universe and beyond

The basic picture of an expanding, cooling Universe leads to a number of startling predictions: the formation of nuclei with its resulting primordial abundances of elements, and the later formation of neutral atoms with the consequent presence of a cosmic background of photons .

As we go back in time, the Universe becomes hotter and hotter and thus the amount of energy available for particle interactions thus increases. As a consequence, the nature of interactions goes from those described at low energy by long-range gravitational and electromagnetic physics, to atomic physics, nuclear physics, all the way to high energy physics at the electroweak scale, gran unification (perhaps), and finally quantum gravity.

The way we know about the high energy interactions of matter is via particle accelerators, which are unravelling the details of these fundamental interactions as we increase their energy. However, one should bear in mind that the physical conditions that take place in our high energy colliders are very different from those that occurred in the early Universe: these machines could never reproduce the conditions of density and pressure in its rapidly expanding thermal plasma. Nevertheless, those experiments are crucial in understanding the nature and rate of the local fundamental interactions available at those energies. What interests cosmologists is the statistical and thermal properties that such a plasma should have, in particular the time at which certain particles decoupled from the plasma, i.e. when their interactions were not quick enough with respect to the expansion of the Universe, and they were left out of equilibrium with the plasma.

One can trace the evolution of the universe from its origin till today. According to the best accepted view, the Universe must have originated at the Planck era ($`10^{19}`$ GeV, $`10^{43}`$ s), from a quantum gravity fluctuation. Needless to say, quantum gravity phenomena are still in the realm of physical speculation, although the last astrophysical probes can already feel effects of quantum gravity at the energy level of the Grand Unified Theories (GUTs), about $`10^{16}`$ GeV ($`10^{35}`$ s)  and beyond . However, it is plausible that a pre Big Bang era of cosmological inflation originated then. The idea of inflation, a period of accelerated expansion before the FRW expansion is established, provides an elegant solution to set up the initial conditions for standard Big Bang cosmology . At some early time, just before the Universe may have reached the GUTs era ($`10^{16}`$ GeV, $`10^{35}`$ s) and was not yet thermalized, the energy density of the Universe was dominated by some material (inflaton field) with negative pressure. Under such an assumption, the Einstein equation leads to the conclusion that during the pre Big Bang era the expansion of the Universe was exponentially accelerated, spreading out the conditions that took place in a small causally connected region over a huge distance across which we observe the Universe now. The last concequence explains, in particular, how became the Universe so big and so uniform before becoming the FRW Universe.

Quantum fluctuations of the inflaton field then left their imprint as tiny perturbations in an otherwise very homogeneous patch of the Universe. One of the most astonishing predictions of inflation is that quantum fluctuations of the inflaton field, being stretched by the exponential expansion generate large-scale perturbations of the space-time geometry. Patterns of perturbations in the geometry are like fingerprints that unequivocally characterize a period of inflation. When matter falls in the gravitational weels of the patterns, it creates density perturbations that collapse gravitationally to form galaxies, clusters and superclusters of galaxies. Perhaps the most intriguing aspect of modern cosmology is the fact that the large-scale structure was seeded at the extremely early epoch of the development of the Universe, when even the conception of matter in usual sense did not exist yet. Moreover this fact is already tested observationally at a very high level of accuracy by the combination of recent results from WMAP satellite  andthe 2dF galaxy redshift survey .

At the end of inflation, the huge energy density of the inflaton field was converted into particles, which soon thermalized and became the origin of the hot Big Bang. Such a process is called reheating of the Universe. Since then, the Universe became radiation-dominated.

It is probable (although by no means certain) that the asymmetry between matter and antimatter originated at the same time as the rest of the energy of the Universe, from the decay of the inflaton. This process is known under the name of baryogenesis , since baryons (mostly quarks at that time) must have originated then, from the leftovers of their annihilation with antibaryons. Any mechanism of baryogenesis requires a violation of the baryon number, C and CP violation, and a departure from thermal equilibrium . The first two conditions can be discussed only within a particle physics model that is definitely beyond the Standard Model. This fact represents perhaps the best example of the perfect interplay between cosmology and particle physics.

As the Universe cooled down, it may have gone through the quark–gluon phase transition ($`10^2`$ MeV, $`10^5`$ s), when baryons (mainly protons and neutrons) formed from their constituent quarks . A minor contribution of antimatter regions, which might be left from baryogenesis, can evolve into condensed antimatter objects; these are important footprints of early phase transitions, which took place far above electroveak energies . In this sense the search for antimatter in space, with the AMS-02 experiment , turns out to be an important step into a deeper understanding of physics beyond the Standard Model which took place in the early Universe.

The furthest window we have on the early Universe at is for the moment that of primordial nucleosynthesis ($`1`$$`0.1`$ MeV, 1 s – 3 min), when protons and neutrons were cold enough that bound systems could form, giving rise to the lightest elements, soon after neutrino decoupling: It is the realm of nuclear physics. The observed relative abundances of light elements are in agreement with the predictions of the hot Big Bang theory. Immediately afterwards, electron–positron annihilation occurs (0.5 MeV, 1 min) and all their energy goes into photons. Much later, at about 1 eV, $`10^5`$ yr, matter and radiation have equal energy densities. Soon after, electrons become bound to nuclei to form atoms (0.3 eV, $`3\times 10^5`$ yr), in a process known as recombination: this is the realm of atomic physics.

Immediately after, photons decouple from the plasma, travelling freely thereafter. Those are the photons we observe as the CMB. The COBE satellite  launched in 1990 observed a tiny anisotropy in the angular distribution of the CMB temperature. It is believed that this anisotropy represents intrinsic fluctuations in the CMB itself, due to the presence of tiny primordial density fluctuations in the cosmological matter present at the time of recombination. These density fluctuations are of great physical interest, because these are the fluctuations that later collapsed to form all of the structures in the Universe, from superclusters to planets. Moreover the angular distribution of the CMB depends, for example, on the baryon density, the Hubble constant $`H_0`$, the densities of dark matter and dark energy. This makes the interpretation of the angular spectrum something of a complex undertaking, but it also makes it a sensitive probe of cosmological models and, implicitly, a particle physics laboratory of high precision .

Much later ($`1`$$`10`$ Gyr), the small inhomogeneities generated during inflation have grown, via gravitational collapse, to become galaxies, clusters of galaxies, and superclusters, characterizing the epoch of structure formation. It is the realm of long–range gravitational physics, dominated by a dark (vacuum) energy. Finally (3 K, 13 Gyr), the Sun, the Earth, and biological life originated from previous generations of stars, and from a primordial soup of organic compounds, respectively.

The recent announcement by the WMAP satellite team of the landmark measurements of the CMB anisotropy  has convincingly confirmed important aspects of the current standard cosmological model described above.

## 4 Conclusions

We have entered a new era in cosmology and astrophysics, where a host of high-precision measurements are already posing challenges to our understanding of the Universe: the density of ordinary matter and the total amount of energy in the Universe; the microwave background anisotropies on a fine-scale resolution; primordial deuterium abundance from quasar absorption lines; the acceleration parameter of the Universe from high-redshift supernovae observations; the rate of expansion from gravitational lensing ; large-scale structure measurements of the distribution of galaxies and their evolution; and many more, which already put constraints on the parameter space of cosmological models.

However, these are only the forerunners of the precision era in cosmology that will dominate the new millennium, and make cosmology a phenomenological science, whose results can be incorporated into experimental high energy physics . An impressive example of such a symbiosis of modern observational cosmology and experimental particle physics was the last improvement of the probe of the neutrino masses by 2dF galaxy redshift surveys data combined with WMAP results. The obtained limit  on the sum ot the neutrino masses of 0.7 eV is substantially better than even the most stringent direct laboratory limit on any individual neutrino mass. On top of it, WMAP data also provide a new limit on the effective number of light neutrino species, beyond the three within the Standard Model : $`1.5<\mathrm{\Delta }N_\nu ^{eff}<4.2`$. This limit is not as stringent as that from LEP, but it applies to additional light particles, that might not be produced in Z decay. The cosmological upper limit on the masses of supersymmetric particles which can play the role of dark matter, is also improved by WMAP .

One of the most difficult challenges that the new cosmology will have to face is understanding the origin of the dark energy, whose behaviour is similar to that of the inflaton field, but 13 Gyr later . It is of tremendous interest from the standpoint of fundamental theory . In this sense cosmology provides us with a way to study a question of central importance for particle theory, namely the nature of the vacuum in quantum field theory. This is something that cannot be studied in particle accelerators, so in this sense cosmology provides a unique window on particle physics.


text:
# Thermal Effects on the Condensates of Chiral Order Parameters, ⟨𝑞̄𝑞⟩ and 𝑔⟨𝑞̄𝜎_{𝜇𝜈}𝐺_{𝜇𝜈}𝑞⟩, and Chiral Restoration from Lattice QCD

## 1 Introduction

The non-perturbative nature of QCD provides various phenomena such as spontaneous chiral-symmetry breaking and color confinement. Currently, the RHIC experiments are in progress to understand these phenomena, and the finite-temperature QCD becomes important. In this point, condensates are useful physical quantities because they characterize the nontrivial QCD vacuum. At finite temperature, the QCD vacuum is expected to change, and we can examine it from the thermal effects on the condensates.

Among various condensates, we study the quark-gluon mixed condensate $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ for the following reasons. First, we note that the mixed condensate is a chiral order parameter independent of $`\overline{q}q`$, since it flips the chirality of the quark as

$`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ (1)
$`=`$ $`g\overline{q}_R(\sigma _{\mu \nu }G_{\mu \nu })q_L+g\overline{q}_L(\sigma _{\mu \nu }G_{\mu \nu })q_R.`$
Therefore, $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$, as well as $`\overline{q}q`$, is an important indicator on the chiral structure of the QCD vacuum. In particular, the thermal effects on both the condensates and the comparison of them will reveal mechanism of the chiral restoration of the QCD vacuum at finite temperature. Second, in contrast to $`\overline{q}q`$, $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ represents a direct correlation between quarks and gluons in the QCD vacuum, and thus characterizes different aspect of the QCD vacuum. Third, the mixed condensate is relevant in various QCD sum rules, especially in baryons , light-heavy mesons  and exotic mesons . In this point, the quantitative estimate of the thermal effects on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ has an impact on hadron phenomenology at finite temperature through the framework of the QCD sum rule. Therefore, it is desirable to estimate $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ at zero/finite temperature by a direct calculation from QCD, such as lattice QCD simulation.

So far, the mixed condensate $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ at zero temperature has been analyzed phenomenologically in the QCD sum rules . In the lattice QCD, a pioneering work  has been done long time ago, but the result was rather preliminary. Recently, new lattice calculations have been performed by us  using the Kogut-Susskind (KS) fermion, and by other group  using the Domain-Wall fermion. At finite temperature, however, there has been no result on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ except for our early results in Ref. . Therefore, we here present the extensive results of the thermal effects on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ as well as $`\overline{q}q`$, including the analysis near the critical temperature.

## 2 Lattice formalism and the parameter set

We calculate the condensates $`\overline{q}q`$ and $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ using the SU(3)c lattice QCD with the KS-fermion at the quenched level. We note that the KS-fermion preserves the explicit chiral symmetry for the quark mass $`m=0`$, which is essential for our study.

The Monte Carlo simulations are performed with the standard Wilson action for $`\beta =6.0,6.1`$ and $`6.2`$. In order to calculate at various temperatures, we use the following lattices as    i) $`\beta =6.0`$$`16^3\times N_t(N_t=16,12,10,8,6,4)`$,  ii) $`\beta =6.1`$$`20^3\times N_t(N_t=20,12,10,8,6)`$,  iii) $`\beta =6.2`$$`24^3\times N_t(N_t=24,16,12,10,8)`$.   We generate 100 gauge configurations for each lattice. However, in the vicinity of the chiral phase transition point, namely, $`20^3\times 8`$ at $`\beta =6.1`$ and $`24^3\times 10`$ at $`\beta =6.2`$, the fluctuations of the condensates get larger. We hence generate 1000 gauge configurations for the above two lattices to improve the estimate. The lattice units are obtained as $`a0.10,0.09,0.07\mathrm{fm}`$ for $`\beta =6.0,6.1,6.2`$, respectively, which reproduce the string tension $`\sigma =0.89\mathrm{GeV}/\mathrm{fm}`$. We calculate the flavor-averaged condensates as

$`a^3\overline{q}q={\displaystyle \frac{1}{4}}{\displaystyle \underset{f}{}}\mathrm{Tr}\left[q^f(x)\overline{q}^f(x)\right],`$ (2)
$`a^5g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ (3)
$`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{f,\mu ,\nu }{}}\mathrm{Tr}\left[q^f(x)\overline{q}^f(x)\sigma _{\mu \nu }G_{\mu \nu }\right],`$
where SU(4)f quark-spinor fields, $`q`$ and $`\overline{q}`$, are converted into spinless Grassmann KS-fields $`\chi `$, $`\overline{\chi }`$ and the gauge-link variable. We adopt the clover-type definition for the gluon field strength $`G_{\mu \nu }`$ on the lattice in order to eliminate $`𝒪(a)`$ discretization error. The more detailed formula for the condensates and gluon field strength are given in Ref. .

In the calculation of the condensates, we use the current-quark mass $`m=21,36,52\mathrm{M}\mathrm{e}\mathrm{V}`$. For the fields $`\chi `$, $`\overline{\chi }`$, the anti-periodic condition is imposed at the boundary in the all directions. We measure the condensates on 16 different physical space-time points of $`x`$ ($`\beta =6.0`$) or 2 points of $`x`$ ($`\beta =6.1,6.2`$) in each configuration. Therefore, at each $`m`$ and temperature, we achieve high statistics in total as 1600 data ($`\beta =6.0`$) or 200 data ($`\beta =6.1,6.2`$). Note that in the vicinity of the critical temperature, $`20^3\times 8`$ ($`\beta =6.1`$) and $`24^3\times 10`$ ($`\beta =6.2`$), we take up to 2000 data so as to achieve high statistics and to guarantee reliability of the results.

## 3 The Lattice results and Discussions

We calculate the condensates $`\overline{q}q`$ and $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ at each quark mass and temperature. We observe that both the condensates show a clear linear behavior against the quark mass $`m`$. Typical plot can be seen in Ref. . Therefore, we fit the data with a linear function and determine the condensates in the chiral limit. Estimated statistical errors are typically less than 5%, while the finite volume artifact is estimated to be about 1%, from the dependence of the condensates on boundary conditions .

For each condensate, We estimate the thermal effects by taking the ratio between the values at finite and zero temperature. This procedure eliminates the renormalization constants because there is no operator mixing for both the condensates in the chiral limit . In figure 1, we plot the thermal effects on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$. We find a drastic change of $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ around the critical temperature $`T_c280\mathrm{M}\mathrm{e}\mathrm{V}`$. This is a first observation of chiral-symmetry restoration in $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$. We also find that the thermal effects on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ are very weak below the critical temperature, namely, $`T\begin{array}{c}<\hfill \\ \hfill \end{array}250\mathrm{M}\mathrm{e}\mathrm{V}`$. The same nontrivial features are also found for $`\overline{q}q`$.

We then compare the thermal effects on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ and $`\overline{q}q`$, because these two condensates characterize different aspects of the QCD vacuum. For this purpose, we plot the thermal effects on $`m_0^2g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q/\overline{q}q`$ in figure 2. From figure 2, we observe that $`m_0^2`$ is almost independent of the temperature, even in the very vicinity of $`T_c`$. This very nontrivial result means that both of $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ and $`\overline{q}q`$ obey the same critical behavior. This might indicate the existence of the universal behavior of chiral order parameters near $`T_c`$, and would impose constraints on the chiral structure of the QCD vacuum.

In summary, we have studied the thermal effects on $`g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q`$ as well as $`\overline{q}q`$ from lattice QCD at the quenched level. We have found a clear signal of chiral restoration from both the condensates, while the thermal effects below $`T_c`$ are very weak. We have also observed that $`m_0^2g\overline{q}\sigma _{\mu \nu }G_{\mu \nu }q/\overline{q}q`$ is almost independent of the temperature, which means both the condensates obey the same critical behavior. This nontrivial similarity would impose constraints on the chiral structure of the QCD vacuum near $`T_c`$. For further studies, a full QCD lattice calculation is in progress in order to analyze the dynamical quark effects on the condensates.


text:
# New results on neutron star low mass transients in the quiescence

## 1 Introduction

Many Low Mass X–ray Binaries (LMXRB) accrete matter at very high rates, and therefore shine as bright X–ray sources only sporadically. Among transient sources are Soft X–ray Transients (SXRTs) host an old neutron star orbiting a low mass star companion (for a review see Campana et al. 1998a). These systems usually alternate periods (weeks to months) of high X–ray luminosity, during which they share the same properties of persistent LMXRBs, to long (1–10 years) intervals of quiescence in which the X–ray luminosity drops by up to 5–6 orders of magnitude.

In recent years the quiescent properties of a handful SXRTs have been studied in some detail thanks to dedicated campaigns carried out with BeppoSAX and RXTE during outburst and observations in quiescence with XMM-Newton and Chandra. Among the main outcomes of these investigations are: 1) SXRT quiescent luminosities are for most cases within a narrow range of $`10^{32}10^{33}\mathrm{erg}\mathrm{s}^1`$ (0.5–10 keV, See Fig. 1); 2) an increasing number of candidate SXRTs in quiescence is being discovered in globular clusters observed with Chandra; 3) the quiescent X–ray spectra of SXRTs display a soft component and in some cases a hard (power-law) component contributing $`1050\%`$ of the flux in the 0.5–10 keV band.

Observationally, the soft component can be modelled with a black body of 0.1–0.3 keV temperature and few km radius. Especially promising is the idea that the soft component of SXRTs may be produced from the cooling of the neutron star heated during repeated outbursts (van Paradijs et al. 1987; Stella et al. 1994; Campana et al 1998a). The theory of deep crustal heating by pycnonuclear reactions compares well with the observations (Brown et al. 1998; Campana et al. 1998a; Rutledge et al. 1999; Colpi et al. 2001). In particular, Rutledge et al. (1999) fitted neutron star atmospheric models to the soft component of quiescent spectra of SXRTs and derived temperatures in the 0.1–0.3 keV and radii consistent with the neutron star radius (10–15 km).

The hard component is well described by a power law tail. In the quiescent spectrum of Aql X-1 and Cen X-4 observed by ASCA and BeppoSAX this component is statistically significant (Asai et al. 1996, 1998; Campana et al. 1998b, 2000) with photon index in the 1–2 range. The same power law is needed (even if less significant statistically) in the analysis of Chandra data in order to achieve an emitting radius of the cooling component consistent with the neutron star radius (otherwise the inferred radius would be smaller; Rutledge et al. 2001a, 2001b). In some sources in globular clusters instead there is no indication off the presence of a hard tail (e.g. Heinke et al. 2003), as well as in a observation of 4U 1608–52 in quiescence with ASCA (Asai et al. 1996). The nature of this hard component is still uncertain. Models range from Comptonization in a hot corona to Advection/Convection Dominate Accretion Flow (ADAF/CDAF; Menou et al. 1999) to emission from a low luminosity jet (Fender et al. 2003) or shock emission from the neutron star that resumed its radio pulsar activity in quiescence. The latter model envisages a situation similar to that of the eclipsing radio pulsar PSR B1259–63 or of the ‘black widow’ pulsar PSR B1957+20: a shock at the boundary between the relativistic MHD wind from the radio pulsar and the matter outflowing from the companion star (Tavani & Arons 1997; Tavani & Brookshaw 1991; Campana et al. 1998a). For the model to explain the observed luminosity in the hard power law component of quiescent SXRTs, a few percent of the pulsar spin-down luminosity must be converted into shock emission. Ongoing deep searches in the radio band have not yet revealed any steady or pulsed emission from quiescent SXRTs (Burgay et al. 2003); however free-free absorption due to matter in the binary system and its surroundings might be an important limiting factor in these searches (Stella et al. 1994; Burgay et al. 2003).

## 2 History

March-April 1997 BeppoSAX observations of Aql X-1 were the first to monitor the evolution of the spectral and time variability properties of a neutron star SXRT from the outburst decay to quiescence. A fast X–ray flux decay was observed, which brought the source luminosity from $`10^{36}`$ to $`10^{33}\mathrm{erg}\mathrm{s}^1`$ in less than 10 days (see Fig. 2). This behaviour was later been observed also in other SXRTs as well as Aql X-1 itself.

The X–ray spectrum showed a pronounced hardening in correspondence of the X–ray flux turn-off (Zhang et al. 1998; Campana et al. 1998b). This result has been confirmed and extended with RXTE data of the 2002 outburst (Campana et al. 2003). At a luminosity of $`10^{35}\mathrm{erg}\mathrm{s}^1`$ the X–ray spectrum showed a power law high energy tail with photon index $`\mathrm{\Gamma }2`$ up to 100 keV, markedly different from the thermal spectrum extending up to 30–40 keV observed at luminosities a factor of 10 higher. At even lower luminosities the X–ray spectrum became even harder with $`\mathrm{\Gamma }11.5`$. These observations, together with the detection by RXTE of a periodicity of a few milliseconds during an X–ray burst, likely indicate that the rapid flux decay is caused by the onset of the propeller effect arising from the very fast rotation of the neutron star magnetosphere (Campana et al. 1998b; Zhang et al. 1998) and that a further transition may occur approaching quiescence (Campana et al. 1998b). This transition may underline the re-activation of a millisecond radio pulsar in quiescence state of Aql X-1, the relativistic wind of which powers the quiescent X–ray emission.

## 3 SAX J1808.4–3658

The study of LMXRBs has been revolutionized by the discovery of the first source (the SXRT SAX J1808.4–3658) showing coherent millisecond pulsations in the persistent emission of its outburst state (Wijnands & van der Klis 1998). Other four millisecond pulsating transients have later been discovered.

SAX J1808.4–3658 showed the same decrease to quiescence as Aql X-1 (Gilfanov et al. 1998) with a hard power law appearing at luminosities of $`5\times 10^{35}\mathrm{erg}\mathrm{s}^1`$ (Gilfanov et al. 1998). A few attempts to reveal SAX J1808.4–3658 in quiescence have been carried out with ASCA and BeppoSAX (Stella et al. 2000; Wijnands et al. 2002a). These observations revealed a source at a luminosity of $`10^{32}\mathrm{erg}\mathrm{s}^1`$; yet the field is crowded, source confusion a problem for the $`1`$ arcmin resolution achieved by BeppoSAX and ASCA and no firm conclusions could be drawn. An XMM-Newton observation with a high signal to noise revealed SAX J1808.4–3658 in quiescence isolated, but a nearby source (about 1 arcmin) can have altered the position and the flux observed with previous instruments (Campana et al. 2002). The 0.5–10 keV unabsorbed luminosity is $`5\times 10^{31}\mathrm{erg}\mathrm{s}^1`$, a relatively low value compared with other neutron star SXRTs. Moreover, at variance with other SXRTs, the quiescent spectrum of SAX J1808.4–3658 was dominated by a hard ($`\mathrm{\Gamma }1.5`$) power law with only a minor contribution ($`\stackrel{<}{}10\%`$), if any, from a soft black body/neutron star atmosphere component (Fig. 3). If the power law originates in the shock between the wind of a turned-on radio pulsar and matter outflowing from the companion, then a spin-down to X–ray luminosity conversion efficiency of $`\eta 10^3`$ is derived; this is in line with the value estimated from the eclipsing radio pulsar PSR J1740–5340 (D’Amico et al. 2001). Within the deep crustal heating model, the faintness of the blackbody-like component indicates that SAX J1808.4–3658 may host a massive neutron star ($`M\stackrel{>}{}1.7M_{}`$). In fact this is the only way to have a rapid cooling thanks to direct Urca process and in turn neutrino cooling.

## 4 Variability in quiescence

Short term variability is potentially a powerful tool for the study of the emission mechanism(s) responsible for the SXRTs quiescent emission. A factor of 3 variability over timescales of days (Campana et al. 1997) and $`40\%`$ over 4.5 yr (Rutledge et al. 2001b) was reported in Cen X-4. Several other neutron star systems have also been found to be variable in quiescence by factors of 3–5 (e.g. Rutledge et al. 2000) but data have been collected over several years and with different instruments. Interesting results have been recently found thanks to the observations of well known systems with Chandra and XMM-Newton.

After years of persistent emission at a high level, KS 1731–260 turned to quiescence (Wijnands et al. 2001; Burderi et al. 2002). The source was observed first with Chandra and 6 months later with XMM-Newton. The quiescent luminosity decrease by a factor of 2–5 (Wijnands et al. 2002b). In the cooling neutron star model, this decrease implies that the crust of the neutron star cooled down rapidly between the two epochs, indicating that the crust has a high conductivity. Further monitoring of KS 1731–260 in quiescence can provide crutial information on the crust conductivity and level of impurities (Rutledge et al. 2002a). A similar turn-off of the X–ray luminosity after years of activity occurred in X1732–304 (in Terzan 1; Wijnands et al. 2002c) and MXB 1659–29 (after 2.5 yr of activity, Wijnands et al. 2002d). These sources may represent a new population of long-lasting transients.

Rutledge et al. (2002b) analyzing Chandra data of the Aql X-1 quiescent phase after the November 2000 outburst found a variable flux and X–ray spectrum. They interpreted these variations in terms of variations of the neutron star effective temperature, which changed from $`130_5^{+3}`$ eV, down to $`113_4^{+3}`$ eV, and finally increased to $`118_4^{+9}`$ eV. Interestingly, during the last observation they also found short-term ($`<10^4`$ s) variability (at $`32\%`$ rms). These data are the first to show an increase of the quiescent flux in the quiescent SXRT. The latter result has a direct impact on the quiescent emission models, since no known mechanism associated with crustal heating can account for this variability. Rutledge et al. (2002b) suggested that accretion might occur during the Aql X-1 quiescent state and variations have can be ascribed to a variable mass inflow rate.

Campana & Stella (2003) reanalyzed the Chandra spectra of Aql X-1, together with a long BeppoSAX observation in the same period, and propose a different interpretation of the spectral variability: that this is due to correlated variations of the power law component and the column density ($`>5`$, a part of which might be intrinsic to the source), while the temperature and flux of the neutron star atmospheric component remained unchanged (Fig. 4). The power law slope vs. column density behaviour is in qualitative agreement with that observed from the radio pulsar binary PSR 1259–63, lending support to the idea that the power law component arises from emission at the shock between a radio pulsar wind and inflowing matter from the companion star.

Thanks to XMM-Newton large throughput, Cen X-4 was observed with the highest signal to noise ever. A $`13`$ hr observation revealed rapid ($`>100`$ s), large ($`45\pm 7\%`$ rms in the $`10^41`$ Hz range) intensity variability, especially at low energies (Fig. 5). In order to investigate the cause of this variability, Campana et al. (2003) divided the data into intensity intervals and fit the resulting spectra with the canonical model for neutron star transients in quiescence, i.e. an absorbed power law plus a neutron star atmosphere model. Variations across different spectra can be mainly accounted for by a variation in the column density together with another spectral parameter (Fig. 6). Based on the available spectra, a variation of the power law could not be preferred over a variation in the temperature of the atmosphere component (even if the first is slightly better in terms of reduced $`\chi ^2`$). This variability can be accounted for by accretion onto the neutron star surface (e.g. Rutledge et al. 2002b) or by the variable interaction between the pulsar relativistic wind and matter outflowing from the companion in a shock front (e.g. Campana & Stella 2003).

## 5 Conclusions

The study of SXRTs in quiescence is revealing a rich phenomenology. For the interpretation of this, different emission mechanism(s) and underlying different physical scenarios are actively debated. These depend on the influence of the neutron star rotation and magnetic field on the infalling matter as well as the physics of deep crustal heating during intense accretion episodes. In the last few years, mainly after the discovery of coherent pulsations in the transient SAX J1808.4–3658, a large interest in these sources has grown. Chandra and XMM-Newton (quiescence) as well as BeppoSAX and RXTE (outburst) observations yielded many interesting results. In particular, the discovery of a sizeable population of candidate SXRTs in globular clusters is providing important new clues on the statistical properties of SXRTs in quiescence. Besides the most studied sources, peculiar behaviours are being found in other sources such as SAX J1808.4–3658 itself, characterized by a factor of 2 lower mean luminosity and, more importantly, a spectrum dominated by the power law component.

On the other hand, in depth studies of as well known sources with Chandra and XMM-Newton are still providing new and unexpected results. In the case of Aql X-1 monitoring its quiescent state with Chandra led to the discovery of luminosity and spectral variability over a month timescale. This can either been interpreted as accretion occurring at these low rates (with an uninfluent magnetic field, Rutledge et al. 2002b) or as correlated variability in the quantity of matter around the system and change in the power law slope due to the interaction of a relativistic pulsar wind with this matter (Campana & Stella 2003). In the case of Cen X-4 an XMM-Newton observation disclosed rapid ($`\stackrel{<}{}100`$ s) X–ray variability, the nature of which, while still uncertain, can be accounted for by the two different mechanisms discussed above (Campana et al. 2003).


text:
# Point-contact study of gap amplitude and symmetry in RuSr2GdCu2O8

## 1 Introduction

The hybrid ruthenate-cuprate superconductor RuSr2GdCu2O8 (Ru-1212) is a triple perovskite material comprising CuO2 bilayers and RuO2 monolayers. It was originally synthesized  with the purpose of incorporating a metallic layer between the CuO2 planes, so as to enhance the inter-layer coupling and reduce the anisotropy. The role of the ruthenium-oxide layers was thus intended to be similar to that of the one-dimensional Cu-O chains in YBCO. In these layers, the ruthenium atoms present the same square-planar coordination as in the CuO2 planes, with a similar bond length. As evidence by X-ray absorption near edge structure (XANES) , nuclear magnetic resonance (NMR)  and magnetization measurements  ruthenium occurs in a mixed valence state as Ru4+ and Ru5+ with concentrations equal to 40% and 60% respectively. As a result, RuO2 planes also act as a charge reservoir, doping the superconducting CuO2 planes.

In addition to the superconducting transition, the compound features magnetic order below the temperature $`T_\mathrm{M}133`$ K, whose nature is still controversial. After the first results suggesting a ferromagnetic order  some theoretical and experimental arguments have been reported in favour of a antiferromagnetic ordering . One possibility to reconcile these opposite evidences is to admit that the ferromagnetically-ordered RuO2 planes are coupled to one another via an antiferromagnetic interaction .

Whatever its nature, the magnetic order is found to survive even when superconductivity sets in. While the coexistence of antiferromagnetism and superconductivity is common among cuprates, the coexistence of superconductivity and ferromagnetism would contrast the widespread belief that these two ordering are competing and mutually exclusive. This is the more true if one takes into account that, according to muon-spin rotation ($`\mu SR`$) measurements , the ferromagnetic order looks uniform and homogeneous, even below $`T_\mathrm{c}`$, on a scale of typically 20Å. This would imply that the interaction between the superconducting and the ferromagnetic order parameter is very weak. As a matter of fact, high-temperature susceptibility data indicate that the ferromagnetic order involves only the Ru moments  while Zn substitutions suggest that only the CuO2 planes host the superconducting charge carriers . The decoupling of CuO2 and RuO2 layers has been recently confirmed by a comprehensive study of DC and AC susceptibility, DC resistance, magnetoresistance , Hall effect and microwave absorption , showing that the RuO2 planes are conducting but do not develop superconductivity.

In spite of the great experimental and theoretical efforts focused on the interplay between superconductivity and magnetism, to our knowledge no spectroscopic studies of the superconducting order parameter in RuSr2GdCu2O8 have been reported so far in literature. In this paper we present and discuss the results of point-contact spectroscopy measurements in polycrystalline RuSr2GdCu2O8 previously characterized by SEM and optical microscopy, AC susceptibility and resistivity measurements. We will show that the order-parameter symmetry that best fits the experimental conductance curves is $`d`$-wave, thus indicating that, as in the other cuprates, superconductivity in Ru-1212 is due to singlet Cooper pairs.

## 2 Experimental details

### 2.1 Preparation and characterization of the samples

The single-phase, polycrystalline samples we used for our measurements were synthesized by solid state reaction starting from high purity stoichiometric powders of RuO2, Gd2O3, CuO and SrCO3. The raw materials were reacted in air at about 960 C, in order to decompose the SrCO3. Then, they were heated in N2 flow at 1010 C and annealed in O2 flow at temperatures ranging from 1050 to 1060 C. Finally, a prolonged anneal in flowing O2 at 1060 C was performed, during which the material densifies, granularity is reduced and ordering within the crystal structure develops. For further details on the preparation technique see Ref..

Figure 1 reports a SEM image of one of the samples studied. It is clearly seen that the thermal treatment  ensures a partial sinterization (some grains begins to coalesce in big aggregates) but a matrix of weakly connected grains is still present. This kind of aggregation can be also perceived by optical inspections with AFM.

Previous to point-contact fabrication and PCS measurements, we quickly characterized the samples from the superconducting and electrical point of view, by means of resistivity and AC susceptibility measurements. Figure 2a reports the temperature dependence of the real and imaginary parts of the susceptibility, measured in the absence of DC magnetic field. The measurements were carried out in a home-made susceptometer designed to allow measurements also in very small samples (sub-mm size). The temperature is given by a SMD-type resistor mounted in the close proximity of the coil containing the sample. Unfortunately, it is calibrated only up to 130 K so that the peak in $`\chi ^{}`$ at about 130 K is only partly visible. The correction for the demagnetizing factor and the correct volume are taken into account in the conversion from the imbalance voltage to the susceptibility. It is clear that, in this kind of measurements, what we actually observe is the volume fraction that is diamagnetically shielded, and *not* the superconducting fraction. In this sense, we can safely say that in our samples the superconducting volume is *at most* 80% of the total volume (see Figure 2a). The diamagnetic transition sets on at $`T_\mathrm{c}^{dia}27`$ K, while the middle point $`T_\mathrm{c}^{mid}`$ (the maximum peak in the imaginary part curves) is estimated at $`T_\mathrm{c}^{mid}23.4`$ K.

The resistivity measurements were carried out by using the conventional four-probe technique and by using AC current to cancel out thermoelectric effects. The voltage signal was recorded by using a EG&G lock-in amplifier that allows using very small frequencies. In particular, the $`\rho (T)`$ curve reported in Figure 2b was obtained with a very small current intensity (3 mA p-p) to reduce possible heating effects, with a frequency of 3.0 Hz. The curve witnesses the good quality of the samples. It is well known that the physical properties of this rutheno-cuprate material – and in particular its critical temperature – are strongly dependent on the details of the preparation procedure. In the present case, the temperature at which the resistance goes to zero ($`T_\mathrm{c}^0=26`$ K) is rather high if compared to other data from literature. The transition sets on at $`T_\mathrm{c}^{onset}=45`$ K, so that $`\mathrm{\Delta }T_\mathrm{c}=19`$ K is the transition width. The change in slope usually attributed to the onset of the magnetic order at $`T_\mathrm{M}`$ is also shown.

### 2.2 Point-contact measurements

Due to the granularity of the material and its only partial sinterization, the surface of the samples looked rather rough and brittle. As a result, point-contact measurements turned out to be rather difficult. We tried with tips of different materials (Au, Pt, Pt-Ir and so on), made by electrochemical etching (with a 25%HNO3\+ 75%HCl solution) starting from wires of different thickness. We also tried to use, instead of the tip, Ag paste or In flakes to make the contacts. In the end, we achieved the best reproducibility of the conductance curves in point contacts obtained by using tips (of either Au or Pt-Ir) made by starting from thick wires ($`\mathrm{}0.1`$ mm), that could apply a rather large pressure on the sample. In these contacts, we were able to follow the temperature evolution of the conductance curves, from 4.2 K up to the critical temperature of the junction.

Once established the contact between tip and sample, we measured the DC $`IV`$ characteristic of the contact and then numerically calculated the differential conductance, i.e. the first derivative d$`I`$/d$`V`$ as a function of the voltage bias $`V`$. To allow the comparison to theoretical models, the resulting curves were then normalized to the normal-state conductance.

Figure 3a reports some examples of normalized conductance curves measured at low temperature ($`T=4.2`$ K). The curves are vertically shifted for clarity. All of them feature a zero-bias conductance peak (ZBCP), and more or less pronounced shoulders at finite voltage. The zero-bias peak is very similar to that produced by Andreev reflection in normal metal/$`d`$-wave superconductor junctions, when the positive interference between incoming and reflected electron gives rise to the zero-energy Andreev bound states. Therefore, the presence of a ZBCP is suggestive of a $`d`$-wave order parameter in Ru-1212.

To investigate in more detail this possibility, we tried to fit the low-temperature normalized conductance curves by using the Blonder-Tinkham-Klapwijk model  generalized to non-conventional order parameter symmetries . For the time being, we only considered the simplest possible symmetries: $`d`$-wave (because of the analogy with other cuprates) and $`p`$-wave (suggested by the comparison with the superconducting ruthenate Sr2RuO4). In Figure 3b an experimental conductance curve (symbols) is reported together with the best-fitting curves obtained with pure $`d`$-wave (solid line) and pure $`p`$-wave (dashed line) order parameter symmetry. The values of the best-fitting parameters are also reported in the legend.

In the case of the $`d`$-wave fit, the free parameters are: the maximum gap amplitude $`\mathrm{\Delta }`$, the barrier transparency parameter $`Z`$ (which is proportional to the potential barrier height at the interface) and the angle $`\alpha `$ the direction of current injection makes with the normal to the interface. In addition to these parameters, a phenomenological broadening parameter $`\mathrm{\Gamma }`$ was needed, that mimics pair-breaking effects due, for example, to non-magnetic impurities. Notice that, however, $`\mathrm{\Gamma }`$ is much smaller than the gap amplitude ($`\mathrm{\Delta }_\mathrm{d}=6.0`$ meV) given by the fit.

As far as the $`p`$-wave symmetry is concerned, we chose one of the triplet $`p`$-wave pair potential symmetries (the so-called $`E_\mathrm{u}`$ states) proposed by Machida *et al*  and by Sigrist and Zhitomirsky . The free fitting parameters are: the gap amplitude $`\mathrm{\Delta }`$, the barrier parameter $`Z`$, the azimuthal angle $`\varphi `$ and, as in the previous case, the broadening parameter $`\mathrm{\Gamma }`$. Interestingly, the amplitude of the $`p`$-wave gap obtained from the fit ($`\mathrm{\Delta }_\mathrm{p}=6.1`$ meV) is consistent with that obtained with the $`d`$-wave symmetry.

It is clearly seen in Figure 3b that the $`d`$-wave fit works much better than the $`p`$-wave one. In particular, the $`d`$-wave curve accounts for both the height and the width of the zero-bias conductance peak, and also follows rather well the shoulders at finite energy. On the contrary, the $`p`$-wave curve has a smoother bell-shaped behaviour that cannot reproduce the main features of the experimental conductance curve. Thus, the $`d`$-wave symmetry is much more compatible with our findings. This fact indicates that the nature of superconductivity in the Ru-1212 compound is by far more similar to other cuprates than to ruthenates. In fact: i) the $`d`$-wave symmetry is common among cuprates, for example in BSCCO; ii) if the order parameter is $`d`$-wave, the Cooper pairs turn out to be spin singlets as in conventional superconductors and HTSC, and not exotic spin-triplet states as in Sr2RuO4. Indirectly, this similarity with other cuprates supports the evidences of a confinement of superconductivity on the CuO2 layers (that are common to Ru-1212 and copper oxides), with no contribution from the Ru-O sublattice .

When the junction was sufficiently stable, we could follow the evolution of its conductance curves on heating, from 4.2 K up to the critical temperature of the junction, where the Andreev-reflection features disappear. Figure 4 reports two examples of temperature dependence of the normalized conductance curves. The curves are vertically shifted for clarity and ordered for increasing values of the temperature. Superimposed to the experimental data (symbols) we also show the best-fitting $`d`$-wave symmetry curves. It is clear that, in general, the fit is satisfactory. In these few cases where the agreement between theoretical and experimental data is poor (see, for example, the upper curve in panel (b)) the choice had to be made between the theoretical curve that best fits the ZBCP and the one that, instead, follows the shoulders.

It is clearly seen that, in both cases, the thermal evolution of the conductance curves is approximately the same, with a progressive decrease of the finite-voltage shoulders and of the zero-bias peak. As a matter of fact, also the behaviour of the best-fit parameters as a function of the temperature is very similar in the two cases. First of all, the broadening parameter $`\mathrm{\Gamma }`$ decreases quickly and goes to zero well before the critical temperature. Second, the gap amplitude starts decreasing at low temperature, contrary to what expected in the conventional BCS theory, and goes to zero at $`T_\mathrm{c}`$ in a sub-linear way. Figure 5 reports the complete $`\mathrm{\Delta }(T)`$ curve for the conductances of Figure 4a, that are better fitted by the model. In the other case, the behaviour is similar but the gap values are slightly greater at intermediate temperatures.

## 3 Conclusions

In conclusion, we have presented some preliminary results of point-contact spectroscopy in the hybrid ruthenate-cuprate superconductor RuSr2GdCu2O8. The low temperature (4.2 K) conductance curves present a clear zero-bias conductance peak, that decreases and goes to zero at the increase of the temperature. The experimental conductance curves turn out to be compatible with a pure $`d`$-wave symmetry of the order parameter, which suggests a tighter similarity of this compound to the other cuprates rather than to the ruthenate Sr2RuO4. It is worth reminding that, however, alternative explanations cannot be ruled out, due to the complexity of the material under study and in particular to its magnetic properties.

One of the authors (V.A.S.) acknowledges the support from RFBR (project N. 02-02-17133) and the Ministry of Science and Technologies of the Russian Federation (contract N. 40.012.1.1.1357).

## References


text:
# Onset and universality of drag reduction in the turbulent Kolmogorov flow

## Abstract

We investigate the phenomenon of drag reduction in a viscoelastic fluid model of dilute polymer solutions. By means of direct numerical simulations of the three-dimensional turbulent Kolmogorov flow we show that drag reduction takes place above a critical Reynolds number $`Re_c`$. An explicit expression for the dependence of $`Re_c`$ on polymer elasticity and diffusivity is derived. The values of the drag coefficient obtained for different fluid parameters collapse onto a universal curve when plotted as a function of the rescaled Reynolds number $`Re/Re_c`$. The analysis of the momentum budget allows to gain some insight on the physics of drag reduction, and suggests the existence of a maximum drag reduction asymptote for this flow.

When a viscous fluid is kept in motion by some external driving, a mean flow is established: the ratio between the work made by the force and the kinetic energy carried by the mean flow is called the drag coefficient, or friction factor. This dimensionless number measures the power that has to be supplied to the fluid to maintain a given throughput. When the flow is laminar, the drag coefficient is inversely proportional to the Reynolds number. Upon increasing the intensity of the applied force the flow eventually becomes turbulent, and the drag coefficient becomes approximately independent of the Reynolds number R1883 , therefore substantially larger than in the viscous case.

In 1949 the British chemist Toms reported that the turbulent drag could be reduced by up to 80% through the addition of minute amounts (few tenths of p.p.m. in weight) of long-chain soluble polymers to water. This observation triggered an enormous experimental activity to characterize this phenomenon (see, e.g., L69 ; V75 ; MC92 ; NH95 ; SW00 ). In spite of these efforts, no fully satisfactory theory of drag reduction is available yet. However, a recent breakthrough has been the observation of drag reduction in numerical simulations of the turbulent channel flow of viscoelastic fluids SBH97 . Most of the features of experimental flows of dilute polymer solutions are successfully reproduced by these models, even at the quantitative level PBNHVH03 . Despite these advances, the understanding of drag reduction in the experimentally relevant geometry of pipe or channel flow is still hindered by the complexity of these flows already at the Newtonian level, i.e. in the absence of polymers DCLPP03 . This consideration motivated us to investigate simpler geometries in the hope that this may shed some light on the basic physical mechanisms of drag reduction (see, e.g., Ref. DCBP02 ).

In this Letter we present the results of an extensive numerical investigation of the viscoelastic turbulent Kolmogorov flow. This system has several analogies with the turbulent channel flow, while its main distinctive trait is the absence of material boundaries. Notwithstanding this major difference we will show that drag reduction takes place in the Kolmogorov flow as well. Furthermore, we observe striking quantitative similarities with experimental results in wall-bounded flows: this points to the conclusion that the basic physical mechanisms of drag reduction be substantially independent of the detailed structure of the flow.

To describe the dynamics of a dilute polymer solution we adopt the linear viscoelastic model (Oldroyd-B) BCAH87

$$_t𝒖+(𝒖\mathbf{})𝒖=\mathbf{}p+\nu _0\mathrm{\Delta }𝒖+\frac{2\eta \nu _0}{\tau }\mathbf{}𝝈+𝑭,$$
(1)
$$_t𝝈+(𝒖\mathbf{})𝝈=(\mathbf{}𝒖)^T𝝈+𝝈(\mathbf{}𝒖)2\frac{𝝈\mathrm{𝟏}}{\tau }+\kappa \mathrm{\Delta }𝝈.$$
(2)
The velocity field $`𝒖`$ is incompressible, the symmetric matrix $`𝝈`$ is the conformation tensor of polymer molecules, and its trace $`\text{tr}𝝈`$ is a measure of their elongation. The parameter $`\tau `$ is the (slowest) polymer relaxation time. The matrix of velocity gradients is defined as $`(\mathbf{}𝒖)_{ij}=_iu_j`$ and $`\mathrm{𝟏}`$ is the unit tensor. The solvent viscosity is denoted by $`\nu _0`$ and $`\eta `$ is the zero-shear contribution of polymers to the total solution viscosity $`\nu =\nu _0(1+\eta )`$. The parameter $`\eta `$ is proportional to the polymer concentration. The diffusive term $`\kappa \mathrm{\Delta }𝝈`$ is added to prevent numerical instabilities SB95 . The constant forcing $`𝑭=(F\mathrm{cos}(z/L),0,0)`$ maintains the system in a statistically stationary state characterized by a mean flow $`𝒖`$. Due to the symmetries of $`𝑭`$, the only nonzero component of the mean velocity is $`u_x`$: it depends on the shear coordinate $`z`$ alone, vanishes at $`z=\pm (\pi /2)L`$, and is even under reflections $`zz`$. Its value at $`z=0`$, $`u_x_{z=0}`$, will be denoted by $`U`$. Finally, we establish a short glossary between the Kolmogorov flow and the channel flow: $`F`$ plays the role of the pressure gradient, $`\pi L`$ is analogous to the channel height, and $`U`$ is equivalent to the centerline velocity.

In this framework, we have performed a series of numerical integrations of eqs. (1) and (2) for a set of values of forcing intensity $`F`$, at fixed $`\nu `$, both for the Newtonian and the viscoelastic case. Comparing results at a given $`F`$ is equivalent to keeping an imposed pressure gradient – therefore a fixed wall-shear stress – in channel flow experiments (see, e.g. Ref. LT88 ). We have measured the mean profiles of several relevant observables, including the average velocity $`u_x`$, the turbulent shear stress (Reynolds stress) $`u_xu_z`$, and the mean polymer stress $`2\nu _0\eta \tau ^1\sigma _{xz}`$. The mean flow is accurately described by the sinusoidal profile $`u_x=U\mathrm{cos}(z/L)`$, both in the Newtonian and in the viscoelastic flow profiles .

However, as shown in Fig. 1, in the latter case the centerline velocity $`U`$ is definitely larger: this is the hallmark of drag reduction. It has to be remarked that – at variance with wall-bounded flows where drag reduction is always accompanied by a structural change in the profile (see e.g. Ref. V75 ) – in the Kolmogorov flow the increase in throughput takes place just by means of an overall rescaling of the mean velocity. This is due to the different boundary conditions: in channel flows, the profile in the viscous sublayer is left unchanged upon polymer addition while the bulk flow increases substantially. This requires a reshaping of the mean profile, that takes actually place through the increase of the extent of the buffer region (see e.g. Ref. LT88 ). In the Kolmogorov flow there is no constraint on velocity profiles, and drag reduction does not necessarily entail their structural change.

To quantify the effect of viscoelasticity on the mean flow, we have defined the drag coefficient as

$$f=\frac{FL}{U^2},$$
(3)
and measured its dependence on the Reynolds number $`Re=UL/\nu `$ Rey-num . For $`Re<\sqrt{2}`$ the flow is laminar with mean velocity $`U=FL^2/\nu `$, giving a drag coefficient $`f=Re^1`$. At $`Re50`$ the system is already in a fully developed turbulent state.

For a Newtonian fluid, numerical data show that the drag coefficient is approximately independent of $`Re`$ (see Fig. 2). This behaviour agrees with the following classical Kolmogorov argument: since the average energy input $`ϵ=FU/2`$ scales as $`ϵ=\frac{\beta }{2}U^3/L`$ in fully developed turbulence, eq. (3) yields a constant drag coefficient $`f=\beta `$. The Newtonian momentum budget gives $`F_x=_zu_xu_z`$ (the viscous contribution being negligible) and therefore a Reynolds stress $`u_xu_z=S\mathrm{sin}(z/L)`$ with $`S=\beta U^2`$. For the turbulent Kolmogorov flow, $`\beta 0.15`$.

When polymers are added $`f`$ may be reduced with respect to its Newtonian value, depending on the polymer elasticity $`El=\nu \tau /L^2`$, the Schmidt number $`Sc=\nu /\kappa `$, and the concentration $`\eta `$, as shown in Fig. 2. For the highest Reynolds number we can attain in our simulations the friction factor is reduced by 75%. Drag reduction is accompanied by changes in the velocity field similar to those occurring in channel flow experiments and simulations: the level of transverse fluctuations $`u_z^2`$ is reduced while longitudinal fluctuations $`(u_xu_x)^2`$ increase and high streamwise velocity streaks are observed (see Fig. 3). Incidentally, we notice that drag reduction is observed at Reynolds numbers definitely smaller than the typical experimental values: this is possible thanks to the relatively high value of elasticity utilized in our simulations. Comparable parameters have been used in numerical simulations of the channel flow as well (see, e.g. Ref. SBH97 ), and produced a similar effect on the threshold for drag reduction.

From the inspection of Fig. 2 we notice some systematic trend: at moderate Reynolds numbers ($`Re200`$) viscoelastic effects do not alter substantially the value of the drag coefficient; at larger $`Re`$ polymers with a higher elasticity are more effective as drag-reducing agents; conversely, polymers with higher diffusivity are less effective. To understand the variation of the drag coefficient with fluid parameters, we sought a dependence of the form $`f=\phi (Re/Re_c)`$ where $`Re_c(El,Sc,\eta )`$ is the critical Reynolds number for the onset of drag reduction. To obtain an explicit expression for $`Re_c`$ we need to extend the argument given in Ref. BFL01  to the case of finite polymer diffusivity. The reasoning goes as follows: for polymers to be substantially elongated, stretching must prevail over elastic relaxation and diffusivity timecr ; at the onset, the terms appearing in eq. (2) must then satisfy $`(u)_c2/\tau +\kappa /L^2`$; since the transition is incipient we can estimate the typical velocity gradient as $`(ϵ_c/\nu )^{1/2}`$, and utilizing the expression $`ϵ_cU_c^3/L`$ we finally obtain

$$Re_c\left(\frac{2}{El}+\frac{1}{Sc}\right)^{2/3}.$$
(4)
For vanishing diffusivity we recover the result of Ref. BFL01 . Extracting the explicit dependence on polymer concentration, we have $`Re_c(1+\eta )^{2/3}`$, which is compatible with the weak dependence on concentration found in experiments V75 .

In Fig. 4 we present the same data as in Fig. 2, now plotted against the rescaled Reynolds number $`Re/Re_c`$. The good quality of the collapse supports the validity of the relation $`f=\phi (Re/Re_c)`$. The function $`\phi `$ is universal with respect to the choice of fluid parameters. Its shape will be derived in the following, with the aid of simple assumptions, starting from the equation for momentum conservation (see Ref. LPPT03  for a similar approach to wall-bounded flows).

Upon time averaging, eq. (1) reduces to $`F_x=\nu _0_z^2u_x+_z(u_xu_z2\nu _0\eta \sigma _{xz}/\tau )`$. Utilizing the numerical observation that the Reynolds stress $`𝒮=u_xu_z=S\mathrm{sin}(z/L)`$ and the polymer stress $`𝒯=2\nu _0\eta \sigma _{xz}/\tau =T\mathrm{sin}(z/L)`$, we obtain the momentum budget $`F=\nu _0U/L^2+S/L+T/L`$. The contribution $`\nu _0U/L^2`$ is relevant only in the laminar regime, and can therefore be neglected. The dependence of the stresses on the rescaled Reynolds number is presented in Fig. 5. Below the threshold the polymer stress is vanishingly small whereas the Reynolds stress is $`S\beta U^2FL`$ in agreement with the observation of a $`Re`$-independent drag coefficient. Above $`Re_c`$, the polymer stress makes a significant contribution to the momentum budget. At the largest $`Re`$ we can attain, the elastic stress reaches almost 50% of the total stress, not far from experimental results PNVH01 . Rescaling the stresses with the critical velocity squared shows that above the onset $`S`$ tends to a constant value $`\gamma U_c^2`$ (see Fig. 5(b)), and the polymer stress follows the law $`T=\delta U^2`$ (Fig. 5(c)). The physical interpretation of these observations is that above the onset of drag reduction an increasing fraction of the momentum injected by the external force is sequestered by polymers, which are however less effective in absorbing it than transverse velocity fluctuations ($`\delta <\beta `$). This results in an enhancement of the mean flow with respect to the Newtonian case, i.e. drag-reduction. Inserting the empirical expressions for $`S`$ and $`T`$, the momentum budget above the onset reads $`F=\gamma U_c^2/L+\delta U^2/L`$, and the resulting drag coefficient is

$$f=\{\begin{array}{ccc}\beta \hfill & \text{for}& ReRe_c,\\ \gamma \left(\frac{Re_c}{Re}\right)^2+\delta \hfill & \text{for}& ReRe_c.\end{array}$$
(5)
This expression is compared with numerical results in Fig. 4, where the values of the parameters $`\gamma `$ and $`\delta `$ have been obtained from the data shown in Fig. 5. The agreement is excellent, except possibly for $`ReRe_c`$, where eq. (5) predicts an abrupt transition: from Fig. 5 this rather appears to be a smooth crossover, whose actual shape cannot be extracted by means of simple arguments. The actual values of $`\beta `$, $`\gamma `$ and $`\delta `$ are not of utmost importance since they are likely to depend on the details of the driving force, and therefore on the shape of the velocity profile. What is crucial to drag reduction is that $`\delta <\beta `$, or – in plain words – that momentum is transferred with greater ease to velocity fluctuations than to elastic ones. Understanding the reasons for this difference would disclose the basic physical mechanisms of drag reduction.

Remarkably, eq. (5) predicts a maximum drag reduction asymptote (see, e.g., Refs. V75 ; SW00 ). Indeed, by increasing the concentration, drag cannot be reduced below the asymptote $`f=\delta `$, independently of polymer elasticity and diffusivity. In this ultimate regime momentum transfer would take place only through polymer stresses. However, the present data do not cover a sufficient span of values of $`Re`$ to allow us to confirm or reject this prediction. Numerical simulations at higher resolution should allow to settle this issue.

We thank M. Chertkov, B. Eckhardt, V. Steinberg, and M. Vergassola for valuable discussions. We acknowledge the support of EU under the contract HPRN-CT-2002-00300 and of MIUR-Cofin 2001023848. Numerical simulations have been performed at CINECA (INFM parallel computing initiative).