format: text

text:
# Andreev reflection off a fluctuating superconductor in the absence of equilibrium

## I Introduction

In the last decades, a considerable effort has been devoted towards the study of the transport properties of normal metal-superconductor (NS) junctions. The situation for the Andreev current and finite frequency noise is well understood when the superconductor is of the BCS type. More recent works have dealt with superconductors whose order parameter has a d-wave symmetry . The role of collective modes arising from the phase of the fluctuations of the superconductor has also been addressed in the framework of linear response theory in high-$`T_c`$ materials , as well as in one-dimensional array of Josephson junctions . Some recent attempts have also included the effect of classical phase fluctuations of the order parameter on Andreev transport, either in the tunneling regime  or using Bogolubov-de Gennes equations .

A physical picture of the Andreev process is drawn in Fig. 1. A voltage bias is applied, such that the chemical potential of the metal $`\mu _N`$ lies $`eV`$ above that of the superconductor $`\mu _S`$. The origin of the energies is taken at $`\mu _S`$. In Andreev reflection with conventional, rigid superconductors, Fig. 1a, an electron below $`\mu _N`$, with an energy $`\mathrm{}\omega `$ is reflected as a hole and a Cooper pair propagates in the superconductor. This pair has zero energy. The reflected hole can be interpreted as an electron below $`\mu _S`$, with energy $`\mathrm{}\omega `$, which enters the superconductor. At zero temperature, the only possible values of $`\mathrm{}\omega `$ are between $`eV`$ and $`eV`$, see Fig. 1a , because both electrons states must have an energy below $`\mu _N`$.

Turning to the case of a superconductor with phase fluctuations, Fig. 1b: exciting a collective mode in the superconductor costs an energy $`\mathrm{}\mathrm{\Omega }`$. One can therefore take an electron of energy $`\mathrm{}(\mathrm{\Omega }/2+\omega )`$ and reflect a hole, which is equivalent to a missing electron of energy$`\mathrm{}(\omega \mathrm{\Omega }/2)`$. The allowed values of $`\mathrm{}\omega `$ must be smaller than $`eV\mathrm{}\mathrm{\Omega }/2`$ (no electrons present above $`\mu _N`$). Note that here, one has to sum over all possible values of $`\mathrm{}\mathrm{\Omega }[0;2eV]`$, taking care of the fact that this is a weighted sum, with a probability distribution $`P(\mathrm{\Omega })`$ which is determined by the Hamiltonian describing the phase.

In this paper, the Andreev current is calculated non-perturbatively in the tunneling Hamiltonian for N-S junctions whose Hamiltonian includes quantum phase fluctuations. The approach is inspired from the Keldysh technique of references , and enables the calculation of the current-voltage characteristics of the Andreev current when the bias voltage $`eV`$ is close to (but smaller than) the gap amplitude $`\mathrm{\Delta }`$. In particular, this corresponds to the experimental situation encountered in scanning tunneling microscope (STM) experiments performed on high-$`T_c`$ materials . Although the calculation is restricted here to an s-wave order parameter, qualitative features can be derived for transport in N-S junctions for arbitrary superconductors, provided that there is a gap. Here, a phase-only effective Hamiltonian will be used, which was previously derived in the literature . In one dimension, this model is isomorphic to an array of Josephson junctions, with Josephson coupling energy $`E_J`$ and charging energy $`E_0`$. The current and the finite frequency noise spectrum $`S(\omega )`$ are both derived, which requires to go beyond linear response theory. In particular, it will be shown that the second derivative of the noise with respect to the frequency $`d^2S/d\omega ^2`$ gives a direct access to information about phase fluctuations, independently of the model chosen to describe the latter.

There are other situations where Andreev scattering is mediated by excitations at the normal supraconductor boundary. If the normal side is replaced by a ferromagnet , the spin of electrons pairs emerging from/entering the superconductor can suffer a spin flip accompanied by the emission/destruction of a magnon. This additional scattering channel effectively may enhance the Andreev current. Except for the case of a half-metal, electrons with opposite spins as well as electrons with the same spin then contribute to the Andreev current. The differences with our phase fluctuations calculation will addressed below.

The paper is organized as follows: in section II , the model Hamiltonian is discussed. In section III, the perturbative scheme which allows to derive the relevant Green’s functions for the transport properties is described. In section IV, the phase fluctuations are introduced in this framework. The current-voltage characteristics associated with the linear Josephson junction array model of Ref.  for the phase appears in section V, as a function of the physical parameters $`E_J`$, $`E_0`$, as well as the barrier transparency.

## II Model Hamiltonian

The Hamiltonian describing the NS junction is composed of three terms: $`H=H_N+H_t+H_{JJA}`$. The first term is simply the Hamiltonian for the metal which is specified to a one band model $`H_N=_{k,\sigma }ϵ_{k,\sigma }c_{k,\sigma }^{}c_{k,\sigma }`$ with width $`W`$. The second term $`H_t`$ is the tunneling term,

$$H_t=\underset{N,S,\sigma }{}\mathrm{\Gamma }c_{N,\sigma }c_{S,\sigma }^{}+h.c.,$$
(1)
where $`\mathrm{\Gamma }`$ is a hopping amplitude which transfers electrons to/from the superconductor.

In the superconductor, it is assumed that the low energy lying excitations are given by fluctuations of the phase of the order parameter. If we forget the poorly screened long-range part of the Coulomb interaction, the Hamiltonian is the same as the one describing an array of Josephson junctions .

$$H_{JJA}=\underset{i,j}{}E_J\mathrm{cos}(\theta _i\theta _j)+\underset{i}{}\frac{q_i^2}{2C},$$
(2)
where $`E_J`$ is the Josephson coupling (between islands $`i`$ and $`j`$), the operator $`\theta _i`$ is the phase of the order parameter at “site” $`i`$, $`q_i`$ is the charge on island $`i`$ and $`\frac{q_i^2}{2C}`$ is the charging energy of island $`i`$. Eq. (2) can be put on a firmer basis, starting from microscopic models. In particular, a phase-only Hamiltonian $`H_S`$ has been derived from various microscopic models such as the attractive Hubbard model :

$$H_S=\underset{i,j}{}\frac{\mathrm{}^2N_S^0}{2m^{}}\mathrm{cos}(\theta _i\theta _j)+\underset{i}{}\frac{1}{2}\frac{m^{}v^2}{N_S^0}\rho _i^2+\underset{i,j}{}\frac{\pi e^2}{\epsilon }\frac{\rho _i\rho _j}{r_ir_j},$$
(3)
where $`N_S^0`$ is the bare superfluid density, $`\epsilon `$ is the dielectric constant, $`r_ir_j`$ is the distance between site $`i`$ and $`j`$. The first term is the Josephson coupling between location $`i`$ and location $`j`$. The second term is the local density fluctuation energy $`(\rho \rho )^2/2\chi _0`$, where $`\rho `$ is the average density and $`\chi _0`$ a density susceptibility. The last term is due to poorly screened long range Coulomb interactions. Since electron and hole-quasiparticles are not quasiparticles for the superconductor, the injected electrons and holes will decay and form pairs that will excite the eigenmodes of $`H_S`$.

## III Non–equilibrium Green’s functions and transport

A voltage bias is introduced between the normal metal and the superconductor and the average current reads

$$I(t)=\frac{ie}{\mathrm{}}\underset{\sigma }{}\mathrm{\Gamma }c_{N,\sigma }^{}(t)c_{S,\sigma }(t)h.c..$$
(4)
Tunneling occurs from the bulk of the normal metal to the bulk of the superconductor through a narrow constriction, the tip of a STM for example.

The Keldysh formalism is used to express the current in terms of Green’s functions. The definitions and properties of this system have been studied in detail in Ref. . Here we adopt the same notations and start from their expression:

$$I(t)=\frac{2e}{\mathrm{}^2}\mathrm{\Gamma }^2_{\mathrm{}}^{\mathrm{}}\left[g_{NN11}^+(t,t_1)G_{SS11}^+(t_1,t)g_{NN11}^+(t,t_1)G_{SS11}^+(t_1,t)\right]𝑑t_1.$$
(5)
The dressed Keldysh Green matrices on the superconducting side are expanded using the Dyson expansion

$$\widehat{G}_{SS}^\pm (t,t^{})=\left[\widehat{G}_{SS}^r\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^\pm \widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^a\right](t,t^{}),$$
(6)
where $`\widehat{\mathrm{\Sigma }}`$ is the self energy matrix describing the hopping from the normal side to the superconductor, and $`\widehat{\mathrm{\Sigma }}^{}`$ describes the opposite process. For simplicity, integrals over time have been omitted. Note that only subgap tunneling is considered here, allowing to set $`\widehat{g}_{SS}^\pm =0`$. Inserting the Dyson expansion, the current reads:

$$\begin{array}{c}I(t)=\frac{2e}{\mathrm{}^4}|\mathrm{\Gamma }|^4_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}G_{SS12}^r(t_1,t_2)G_{SS21}^a(t_3,t)\times \hfill \\ \times \left(g_{NN11}^+(t,t_1)g_{NN22}^+(t_2,t_3)g_{NN11}^+(t,t_1)g_{NN22}^+(t_2,t_3)\right)dt_1dt_2dt_3.\hfill \end{array}$$
(7)
The dressed, retarded Green’s function on the superconducting side is given by:

$$G_{SS12}^r(t^{},t)=g_{SS12}^r(t^{},t)+_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\left(\widehat{g}_{SS}^r(t^{},t_1)\widehat{\mathrm{\Sigma }}_S^r(t_1,t_2)\widehat{G}_{SS}^r(t_2,t)\right)_{12}𝑑t_1𝑑t_2+\mathrm{},$$
(8)
with the dressed self energy $`\widehat{\mathrm{\Sigma }}_S^r`$ defined as:

$$\widehat{\mathrm{\Sigma }}_S^r(t_1,t_2)=_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\widehat{\mathrm{\Sigma }}^{}(t_1,t_1^{})\widehat{g}_{NN}^r(t_1^{},t_2^{})\widehat{\mathrm{\Sigma }}(t_2^{},t_2)𝑑t_1^{}𝑑t_2^{},$$
(9)
and the unperturbed momentum averaged Green’s function in the superconductor reads :

$$g_{SS12}^r(\omega )=\frac{\mathrm{}}{W}\frac{\mathrm{\Delta }}{\sqrt{\mathrm{\Delta }^2(\mathrm{}\omega +i\alpha )^2}}e^{i\theta _0}\stackrel{~}{g}_{SS12}^r(\omega )e^{i\theta _0},$$
(10)
where $`\alpha `$ is an infinitesimally small positive quantity and $`\theta _0`$ the phase of the superconductor, which appears also in the gap parameter: $`\mathrm{\Delta }=\mathrm{\Delta }e^{i\theta _0}`$.

## IV Transport in the presence of phase fluctuations

In the presence of phase fluctuations, the Green’s function on the superconducting side contains two time scales . First, fast time scales (of the order $`\mathrm{}/\mathrm{\Delta }`$) which were described in the previous section in Eq. (10). Secondly, “slow” time scales associated with the fluctuating phase $`\theta (x,t)`$:

$$\mathrm{\Delta }=\mathrm{\Delta }e^{i\theta (x,t)}.$$
(11)
The phase of the order parameter has both spatial and temporal variations and is no longer equal to $`\theta _0=\theta (0,0)`$. Tunneling is supposed to occur close to the interface. The typical dimension parallel to the interface will be small with respect to both the classical correlation length $`\xi _\varphi ^c`$  and to the amplitude correlation length $`\xi _a`$. In high-$`T_c`$ superconductors, $`\xi _\varphi ^c`$ and $`\xi _a`$ are of the order of a few nanometers. $`\xi _\varphi ^c`$ is believed to be generally not much larger than $`\xi _a`$, and $`\xi _a`$ is much smaller than in BCS superconductors. We shall thus replace $`x`$ by zero. The temporal variations of $`\theta `$ are due to quantum fluctuations; the phase of the order parameter and the number of pairs are conjugated variables. These variations occur on time scales of the order at least of the plasma frequency $`\omega _p\sqrt{8E_0E_J}/\mathrm{}`$ with $`E_0=q_0^2/2C`$ and $`E_J`$ appear in the definition of the Hamiltonian of Eq. (2). Equivalently, $`\mathrm{}\omega _p<\mathrm{\Delta }`$, which is the contrary to what usually occurs in BCS superconductors. The dynamic nature of the phase is reflected in the anomalous Green’s functions :

$$g_{SS12}^{r,a}(t^{},t)=\stackrel{~}{g}_{SS12}^{r,a}(t^{},t)e^{i\frac{1}{2}\left(\theta (0,t)+\theta (0,t^{})\right)}$$
(12)
where $`\stackrel{~}{g}_{SS12}^{r,a}(t^{},t)`$ is the inverse Fourier transform of $`\stackrel{~}{g}_{SS12}^{r,a}(\omega )`$, and $`g_{SS21}^{r,a}(t^{},t)`$ has the opposite phase. Because of the slow temporal variation of the phase, the exponential in Eq. (12): $`\mathrm{exp}[i(\theta (0,t)+\theta (0,t^{}))/2]\mathrm{exp}[i\theta (0,t^{})]`$ without affecting the results in a significant manner.

### A Current and differential conductance

Substituting Eq. (12) in the current, the fourth (lowest non vanishing) order term in $`\mathrm{\Gamma }`$ which contributes to the current is:

$$\begin{array}{c}I(t)=\frac{2e}{\mathrm{}^4}\mathrm{\Gamma }^4_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\stackrel{~}{g}_{SS21}^a(t_3,t)\stackrel{~}{g}_{SS12}^r(t_1,t_2)e^{i\theta (0,t_1)}e^{i\theta (0,t_3)}\hfill \\ \times \left(g_{NN11}^+(t,t_1)g_{NN22}^+(t_2,t_3)g_{NN11}^+(t,t_1)g_{NN22}^+(t_2,t_3)\right)dt_1dt_2dt_3.\hfill \end{array}$$
(13)
The bracket $`\mathrm{}`$ signifies that an average has been performed over the dynamical degrees of freedom of the phase. Depending on the nature of the phase, this will correspond to a Gaussian model – the so called ordered phase – or to a nonlinear model  (see below). Note that this expression is analogous to the one derived several decades ago by Combescot et al. for inelastic tunneling : there the working assumption was that the inelastic coupling did not couple directly the two electrodes. Here, this same assumption is natural because the phase fluctuations happen only on the superconducting side. The correlator of the phase and its Fourier transform are denoted:

$`p(t)`$ $`=`$ $`e^{i\theta (0,t)}e^{i\theta (0,0)},`$ (15)
$`P(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{i\omega t}p(t)𝑑t.`$ (16)
$`P(\omega )`$ can also be viewed as the occupation probability of the collective phase modes with energy $`\mathrm{}\omega `$. Using the Fourier representation, the current becomes:

$$\begin{array}{c}I(t)=\frac{2e}{h\mathrm{}^3}\mathrm{\Gamma }^4_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\stackrel{~}{g}_{SS21}^a(\omega +\mathrm{\Omega }/2)\stackrel{~}{g}_{SS12}^r(\omega \mathrm{\Omega }/2)\hfill \\ \times \left[g_{NN11}^+(\omega +\mathrm{\Omega }/2)g_{NN22}^+(\omega \mathrm{\Omega }/2)g_{NN11}^+(\omega +\mathrm{\Omega }/2)g_{NN22}^+(\omega \mathrm{\Omega }/2)\right]P(\mathrm{\Omega })d\mathrm{\Omega }d\omega .\hfill \end{array}$$
(17)
Next, for frequencies below the gap, one can approximate $`\stackrel{~}{g}_{SSij}^{a,r}(\omega \pm \mathrm{\Omega }/2)\stackrel{~}{g}_{SSij}^{a,r}(\omega )\mathrm{}/\left(W\sqrt{1\mathrm{}^2\omega ^2/|\mathrm{\Delta }|^2}\right)`$ ($`W`$ is the bandwidth), and the current becomes:

$$\begin{array}{c}I(t)=\frac{2e}{h\mathrm{}^3}\mathrm{\Gamma }^4_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\stackrel{~}{g}_{SS21}^a(\omega )\stackrel{~}{g}_{SS12}^r(\omega )\hfill \\ \times \left[g_{NN11}^+(\omega +\mathrm{\Omega }/2)g_{NN22}^+(\omega \mathrm{\Omega }/2)g_{NN11}^+(\omega +\mathrm{\Omega }/2)g_{NN22}^+(\omega \mathrm{\Omega }/2)\right]P(\mathrm{\Omega })d\mathrm{\Omega }d\omega ,\hfill \end{array}$$
(18)
where the two causal Green function elements in Eq. (18) are related to the electron density of states $`\rho (\mathrm{}\omega )`$:

$`g_{NN11}^+(\omega )`$ $`=`$ $`2\pi i\widehat{\rho }_{11}(\mathrm{}\omega eV)n_F(\mathrm{}\omega eV),`$ (20)
$`g_{NN22}^+(\omega )`$ $`=`$ $`2\pi i\widehat{\rho }_{22}(\mathrm{}\omega +eV)\left[n_F(\mathrm{}\omega +eV)1\right],`$ (21)
with $`\widehat{\rho }_{11}(\mathrm{}\omega eV)\widehat{\rho }_{22}(\mathrm{}\omega +eV)\rho (\mathrm{}\omega )\mathrm{}/\pi W`$.

Using the property $`\stackrel{~}{g}_{SS21}^a(\omega )=[\stackrel{~}{g}_{SS12}^r(\omega )]^{}`$, the current in the absence of phase fluctuation is recovered by setting $`P(\mathrm{\Omega })`$ to be a delta function:

$$I_0=\frac{8e}{h\mathrm{}}\frac{\mathrm{\Gamma }^4}{W^2}_{\mathrm{}}^{\mathrm{}}\stackrel{~}{g}_{SS21}^a(\omega )^2\left[n_F(\mathrm{}\omega eV)n_F(\mathrm{}\omega +eV)\right]𝑑\omega ,$$
(22)
where subscript $`0`$ means that in this limit $`E_0/E_J=0`$. The known result for the Andreev current is recovered . It is well known however that even in the absence of phase fluctuations, a perturbative calculation to fourth order in $`\mathrm{\Gamma }`$ is not sufficient to give a satisfactory answer for the whole range of voltage biases between $`0`$ and $`\mathrm{\Delta }`$. Indeed, a resummation procedure, explicited by Eq. (8), has to be carried out. In the absence of phase fluctuations, it is thus sufficient to replace $`\stackrel{~}{g}_{SS21}^r(\omega )`$ by $`G_{SS21}^r(\omega )`$: this allows to recover the scattering theory results of BTK .

Turning back to the lowest order contribution of the current in the presence of phase fluctuations, one obtains:

$$\begin{array}{c}I=\frac{8e}{h\mathrm{}}\frac{\mathrm{\Gamma }^4}{W^2}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}G_{SS21}^a(\omega )G_{SS12}^r(\omega )\hfill \\ \times \left[n_F(\mathrm{}\omega eV+\mathrm{}\mathrm{\Omega }/2)n_F(\mathrm{}\omega +eV\mathrm{}\mathrm{\Omega }/2)\right]P(\mathrm{\Omega })d\mathrm{\Omega }d\omega .\hfill \end{array}$$
(23)
Note that $`P(\mathrm{\Omega })`$ decreases when $`\mathrm{\Omega }`$ becomes large. In particular, the integrated contribution $`_\mathrm{\Omega }^{\mathrm{}}P(\mathrm{\Omega }^{})𝑑\mathrm{\Omega }^{}`$ becomes negligible well before $`\mathrm{}\mathrm{\Omega }>2eV`$, allowing modify the upper bound of the integral over $`\mathrm{\Omega }`$. Simultaneously, the advanced (retarded) Green’s functions are dressed as before and take the form :

$$G_{SS12}^{r,a}(\omega )=\frac{1}{D^{r,a}(\omega )}\stackrel{~}{g}_{SS12}^{r,a}(\omega ),$$
(24)
so that at zero temperature the current reads:

$$I\frac{16e|\mathrm{\Gamma }|^4}{h\mathrm{}W^4}\mathrm{\Delta }^2_0^{2\frac{eV}{\mathrm{}}}P(\mathrm{\Omega })𝑑\mathrm{\Omega }_0^{\frac{eV}{\mathrm{}}\frac{\mathrm{\Omega }}{2}}\frac{1}{D(\omega )}𝑑\omega ,$$
(25)
with

$`D(\omega )`$ $`=`$ $`\left({\displaystyle \frac{|\mathrm{\Delta }|^2}{\mathrm{}^2}}\omega ^2\right)D^r(\omega )D^a(\omega )=\left({\displaystyle \frac{|\mathrm{\Delta }|^2}{\mathrm{}^2}}\omega ^2\right)|D^a(\omega )|^2`$ (26)
$`=`$ $`\left({\displaystyle \frac{\mathrm{\Delta }^2}{\mathrm{}^2}}\omega ^2\right)\left(1+{\displaystyle \frac{|\mathrm{\Gamma }|^4}{W^4}}\right)^2+4{\displaystyle \frac{|\mathrm{\Gamma }|^4}{W^4}}\omega ^2.`$ (27)
In a more rigorous approach, the expansions of $`G_{SS12}^r(\omega )`$ have to be performed in the presence of phase fluctuations. This means that they should contain not only linear terms in $`P(\mathrm{\Omega })`$ but also higher order correlators of the exponentiated phases. Here, in the expression for the current in Eq. (25), only the first order in $`P(\mathrm{\Omega })`$ has been retained, which constitutes the analog of a weak inelastic coupling – single phonon – approximation of Ref. .

The differential conductance $`dI/dV`$ can be straightforwardly computed from Eq. (25):

$$\frac{dI}{dV}\frac{16e^2|\mathrm{\Gamma }|^4}{h\mathrm{}^2W^4}\mathrm{\Delta }^2_0^{2\frac{eV}{\mathrm{}}}P(y)\frac{1}{D\left(\frac{eV}{\mathrm{}}\frac{y}{2}\right)}𝑑y.$$
(28)
### B Finite frequency noise

Another useful tool which probes the fluctuations is the finite frequency noise. The symmetrized real time noise correlator is given by:

$$S(tt^{})=\delta I(t)\delta I(t^{})+\delta I(t^{})\delta I(t)$$
(29)
where $`\delta I(t)=I(t)I(t)`$. This current–current correlator is expressed in terms of the Keldysh Green’s functions: using the self energy matrix, we can write it as a trace.

$`S(t,0)`$ $`=`$ $`e^2\mathrm{Tr}[\{(\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SN}^+(0,t)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SN}^+(t,0)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}^{}\widehat{G}_{NN}^+(0,t)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^+(t,0))+(\widehat{\mathrm{\Sigma }}\widehat{\mathrm{\Sigma }}^{};SN)\}`$ (31)
$`+\{(0,t)(t,0)\}],`$
where $`\widehat{\sigma }_z`$ is a Pauli matrix.

The lowest non vanishing contribution (order 4 in $`\mathrm{\Gamma }`$) is extracted by expanding the Green’s functions in Eq. (31). The corresponding Dyson equations read (time integrations are implicit):

$`\widehat{G}_{SN}^\pm `$ $`=`$ $`\widehat{G}_{SS}^r\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^\pm ,`$ (32)
$`\widehat{G}_{NS}^\pm `$ $`=`$ $`\widehat{g}_{NN}^\pm \widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^a,`$ (33)
$`\widehat{G}_{NN}^\pm `$ $`=`$ $`\widehat{g}_{NN}^\pm ,`$ (34)
where we have purposely written the Green’s functions to the lowest order which they contribute. Inserting these expressions in the real time correlator gives:

$`S(t,0)=e^2{\displaystyle _{\mathrm{}}^+\mathrm{}}dt_1{\displaystyle _{\mathrm{}}^+\mathrm{}}dt_2\mathrm{Tr}\left[\right\{`$ $`\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^r(0,t_1)\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(t_1,t)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^r(t,t_2)\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(t_2,0)`$ (35)
$`+`$ $`\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(0,t_1)\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^a(t_1,t)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(t,t_2)\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^a(t_2,0)`$ (36)
$``$ $`\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^r(0,t_1)\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(t_1,t_2)\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^a(t_2,t)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(t,0)`$ (37)
$``$ $`\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(0,t)\widehat{\sigma }_z\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^r(t,t_1)\widehat{\mathrm{\Sigma }}^{}\widehat{g}_{NN}^+(t_1,t_2)\widehat{\mathrm{\Sigma }}\widehat{G}_{SS}^a(t_2,0)\}+\{0t\}].`$ (38)
The procedure for including the phase fluctuations in the noise correlator is identical to that used in the expression of the current of Eq. (7). However, here additional contributions which are proportional to $`\widehat{G}_{SS}^r\widehat{G}_{SS}^r`$ or to $`\widehat{G}_{SS}^a\widehat{G}_{SS}^a`$ occur, but these turn out not contribute to the noise to lowest order in the phase correlator.

The current noise spectrum $`S(\omega )`$ is the Fourier transform of $`S(t)`$. In the subgap regime, the only non–vanishing elements of the causal Green’s functions are the off diagonal ones : $`\stackrel{~}{g}_{SS12}^{a,r}=\stackrel{~}{g}_{SS21}^{a,r}\mathrm{}/W`$. The current noise spectrum can be expressed in terms of the off diagonal Keldysh Green’s functions:

$`S(\omega )`$ $``$ $`{\displaystyle \frac{4e^2}{\mathrm{}^2}}{\displaystyle \frac{|\mathrm{\Gamma }|^4}{2\pi W^2}}{\displaystyle _0^{2\frac{eV}{\mathrm{}}}}𝑑\mathrm{\Omega }P(\mathrm{\Omega }){\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega _1`$ (43)
$`[g_{N22}^+(\omega _1\mathrm{\Omega })g_{N11}^+(\omega _1\omega )`$
$`+g_{N11}^+(\omega _1+\mathrm{\Omega })g_{N22}^+(\omega _1\omega )`$
$`+g_{N11}^+(\omega _1+\mathrm{\Omega })g_{N22}^+(\omega _1+\omega )`$
$`+g_{N22}^+(\omega _1\mathrm{\Omega })g_{N11}^+(\omega _1+\omega )].`$
The integration over $`\omega _1`$ defines the relevant energy intervals which contribute to the noise:

$`S(\omega )`$ $``$ $`{\displaystyle \frac{16e^2}{2\pi }}{\displaystyle \frac{|\mathrm{\Gamma }|^4}{W^4}}\left[{\displaystyle _0^{2\frac{eV}{\mathrm{}}\omega }}𝑑\mathrm{\Omega }P(\mathrm{\Omega })\mathrm{\Theta }\left({\displaystyle \frac{2eV}{\mathrm{}}}\omega \mathrm{\Omega }\right)\left({\displaystyle \frac{2eV}{\mathrm{}}}\omega \mathrm{\Omega }\right)+{\displaystyle _0^{2\frac{eV}{\mathrm{}}+\omega }}𝑑\mathrm{\Omega }P(\mathrm{\Omega })\left({\displaystyle \frac{2eV}{\mathrm{}}}+\omega \mathrm{\Omega }\right)\right],`$ (44)
where $`\mathrm{\Theta }(\mathrm{\Omega })`$ is the Heaviside function.

In the limit of low frequencies, the noise becomes:

$`S(\omega =0)`$ $`=`$ $`{\displaystyle \frac{32e^2}{2\pi }}{\displaystyle \frac{|\mathrm{\Gamma }|^4}{W^4}}{\displaystyle _0^{2\frac{eV}{\mathrm{}}}}𝑑\mathrm{\Omega }P(\mathrm{\Omega })\left({\displaystyle \frac{2eV}{\mathrm{}}}\mathrm{\Omega }\right)`$ (46)
$`=`$ $`4eI,`$ (47)
which corresponds to the Schottky formula. Note that this constitutes an illustration of Schottky formula for a situation where the charge transfer is effectively inelastic: as shown in Fig. 1b, the transfer of two electrons in the superconductor generates phase quanta and thus constitutes an inelastic process. These processes are included when the distribution $`P(\omega )`$ in the current and noise deviates from a delta function.

The current noise spectrum can be measured experimentally and provides a direct way to obtain information on $`P(\omega )`$ and thus on phase fluctuations. In the tunnel limit, the distribution $`P(\omega )`$ is related to the second derivative of the current noise spectrum by the formula

$$\frac{d^2S}{d\omega ^2}\frac{\mathrm{}I_0}{V}\left[P\left(\frac{2eV}{\mathrm{}}+\omega \right)+P\left(\frac{2eV}{\mathrm{}}\omega \right)\mathrm{\Theta }\left(\frac{2eV}{\mathrm{}}\omega \right)\right].$$
(48)
The prefactor stands for the current in the absence of phase fluctuations :$`I_0=(16e^2V|\mathrm{\Gamma }|^4)/(\mathrm{}W^4)`$.

Note that so far, no specific model for the phase fluctuations has been specified, as the only assumption used what the fact that $`P(\mathrm{\Omega })`$ falls sufficiently rapidly under the gap. In the next section, current voltage characteristics are plotted using a specific model to describe the phase dynamics.

## V Application to the 1D Josephson-Junction Array model

One spatial dimension is specified, motivated by the fact that tunneling typically occurs near the tip of a scanning tunneling microscope. The dynamics of the phase can then be modeled by a one-dimensional Josephson-junction Array (JJA), which was previously studied in linear response . The phase correlator in Eq. (15) is then known. The parameters of this model are the Josephson energy $`E_J`$ and the charging energy $`E_0`$ (which in turn define the plasma frequency), together with the parameter $`\kappa `$ which depends on the ratio of the junction capacitance divided by the capacitance of the islands. If the latter capacitance dominates over that of the junction, $`\kappa =2`$.

The importance of phase fluctuations is monitored by the ratio $`\mathrm{}\omega _p/E_J`$. In particular, for $`\mathrm{}\omega _p/E_J<\pi /2`$, the phase Hamiltonian can be mapped to a Luttinger liquid model and $`P(\mathrm{\Omega })`$ can be derived in a standard way from harmonic oscillator correlators. On the other hand, for $`\mathrm{}\omega _p/E_J>\pi /2`$, the phase correlator decays exponentially in time. This is identified as the disordered phase. At $`\mathrm{}\omega _p/E_J=\pi /2`$, there is a Kosterlitz-Thouless transition between the ordered and the disordered phase.

### A Ordered phase

At zero temperature, the time dependent phase correlator is given by:

$$p(t)=\mathrm{exp}\left\{\frac{\sqrt{2}}{\pi }\sqrt{\frac{E_0}{E_J}}\left[_0^{\kappa \omega _pt}\frac{1\mathrm{cos}x}{x}\frac{dx}{\sqrt{1(x/\kappa \omega _pt)^2}}+i\mathrm{sgn}(t)_0^{\kappa \omega _pt}\frac{\mathrm{sin}x}{x}\frac{dx}{\sqrt{1(x/\kappa \omega _pt)^2}}\right]\right\},$$
(49)
with $`\mathrm{sgn}(t)`$ the sign function. Its Fourier transform $`P(\mathrm{\Omega })`$ is identically zero for $`\mathrm{\Omega }<0`$. For $`\mathrm{\Omega }>0`$, it has a power law behavior close to $`\mathrm{\Omega }=0`$:

$$P(\mathrm{\Omega })\mathrm{\Omega }^{\left(1\frac{1}{\pi }\sqrt{\frac{2E_0}{E_J}}\right)},$$
(50)
and goes rapidly to zero for $`\mathrm{\Omega }>\kappa \omega _p`$.

For $`\mathrm{\Omega }=\kappa \omega _p`$, $`P(\mathrm{\Omega })`$ shows a weaker singularity  and diverges as $`(\omega _p\mathrm{\Omega })^\alpha `$, with $`\alpha =1/2\pi ^1\sqrt{2E_0/E_J}`$ (here this corresponds to fixing the parameter of Ref.  $`\kappa =2`$).

Turning to the transport properties, the effect of the phase fluctuations on the noise are probed. For convenience, we compare our results to those derived with a BTK model, which corresponds to an interface with a rigid superconductor. The quantity $`\gamma =\mathrm{\Gamma }/W`$ corresponds to the transparency of the barrier. $`\gamma `$ enters the prefactor of the delta function potential of BTK theory with the dependence $`\mathrm{}^2k_F(1\gamma ^2)/2m\gamma `$ ($`k_F`$ is the Fermi wave vector). $`\gamma =1`$ corresponds to a perfect contact, while $`\gamma 0`$ describes an opaque barrier.

In Fig. 2 the second derivative of the symmetrized noise spectrum $`d^2S/d\omega ^2`$ is plotted as a function of $`\omega `$, for a perfectly good contact. Several values of the parameter $`\mathrm{}\omega _p/E_J`$ are considered, starting from weak fluctuations to substantial collective excitations within the ordered phase. To illustrate the effect of phase fluctuations, the choice $`\mathrm{}\omega _p=|\mathrm{\Delta }|/2`$ is made, so that all the features in $`P(\mathrm{\Omega })`$ are located within the superconducting gap.

In the absence of phase fluctuations, the finite frequency noise has a singular derivative at $`\mathrm{}\omega =2eV`$ . In Fig. 2, this would imply a delta function peak in $`d^2S/d\omega ^2`$. For low voltages, there is no deviation with respect to the results of scattering theory. The effect of phase fluctuations is twofold: first, the delta function broadens to a power law divergence (partially cut on the figure), and second, it acquires a secondary peak – previously discussed for $`P(\mathrm{\Omega })`$ – at $`2eV\kappa \mathrm{}\omega _p`$. This illustrates the result of Eq. (48): an experimental measurement of the current noise spectrum constitutes a direct diagnosis of the importance of phase fluctuations.

Naturally, the differential conductance $`dI/dV`$ in Eq. (28) is also affected by the phase fluctuations. It is studied below for different regimes: first, it is plotted as a function of bias, for both the case of a perfect contact and for the case of a weak transmission; second, the transparency is varied while the bias is fixed to a large value in order to explore the deviations from the tunnel limit. In each case, curves are obtained for several values of the ratio $`\mathrm{}\omega _p/E_J`$ which characterizes the importance of the fluctuations.

The differential conductance of a perfect junction is displayed in Fig. 3 as a horizontal line at $`dI/dV=4e^2/h`$ . In the presence of fluctuations, it saturates to the BTK value for biases larger than $`\mathrm{}\omega _p/2`$ (here the bias is varied from $`0`$ to $`|\mathrm{\Delta }|`$): for sufficiently large $`eV`$, $`_0^{2eV/\mathrm{}}P(\mathrm{\Omega })𝑑\mathrm{\Omega }`$ is essentially equal to unity because for $`\gamma =1`$, $`D(x)`$ which enters Eq. (28) is a constant. Phase fluctuations have a tendency to decrease the differential conductance at low bias. In particular for $`\mathrm{}\omega _p/E_J1`$, the significant contributions from the integral $`_{\omega _p}^{2eV/\mathrm{}}P(\mathrm{\Omega })𝑑\mathrm{\Omega }`$ are the cause for the deviations from the ideal conductance at large bias.

The case of $`\gamma =1`$ is somewhat academic because perfect contacts are difficult to achieve in practice. The same curves are plotted now for $`\gamma =0.4`$ in Fig. 4. For intermediate biases, phase fluctuations do not affect the differential conductance drastically: at $`eV=\mathrm{}\omega _p`$ for instance, $`dI/dV`$ is essentially constant (of the order $`4e^2\gamma ^2/h`$), and all curves can not be distinguished. However, deviations occur both at small and large voltages. For small voltages, according to Eq. (28) and due to the fact that $`P(0)=0`$ the differential conductance in the presence of any phase fluctuation is required to vanish (this is not really obvious in Fig. 4 because of the choice of scale). More dramatic is the fact that phase fluctuations play a role at large voltages, contrary to the high transparency regime.

In Fig. 5, the transparent regime is examined with a bias voltage $`eV=|\mathrm{\Delta }|`$. In contrast to the previous curves, the differential conductance is plotted versus the transparency of the barrier. The plasma frequency and the ratio $`\mathrm{}\omega _p/E_J`$ are chosen as previously.

For barrier transparencies smaller than $`0.1`$, the differential conductance has power law dependence ($`dI/dV\gamma ^4`$). As the transparency is increased, in the presence of phase fluctuations an inflexion point appears and the rise of $`dI/dV`$ is then slower. Note once again that it is reduced from the BTK case when phase fluctuations become important. As indicated previously in Fig. 3, the differential conductance saturates to the BTK value at high $`\gamma `$ regardless of the degree of phase fluctuations.

### B Disordered phase

In the disordered phase ($`\mathrm{}\omega _p/E_J>\pi /2`$), the behavior of the Fourier transform of the time dependent phase correlator close to the transition has been previously derived from the $`XY`$ model . For small $`\mathrm{\Omega }`$,

$$P(\mathrm{\Omega })\mathrm{\Theta }(\mathrm{\Omega }\omega _p\xi ^1)(\mathrm{\Omega }\omega _p\xi ^1)^{1/2},$$
(51)
where $`\mathrm{\Theta }(x)`$ is the Heaviside function. The results of linear response theory of Ref.  are briefly recalled and compared to the predictions of the Keldysh calculation of Eq. (25). In this phase, a threshold voltage $`\mathrm{}\omega _p\xi ^1`$ has to be reached in order to obtain a non zero current. Here, the Kosterlitz-Thouless correlation length (in dimensionless units) reads $`\xi =\mathrm{exp}\left[b\left(2/\pi \sqrt{E_J/8E_0}\right)^{1/2}\right]`$, where $`b`$ is a positive constant. Because the goal is to predict the current as a function of voltage for both the tunnel limit and the transparent regime (for $`0<eV<\mathrm{\Delta }`$), $`P(\mathrm{\Omega })`$ needs to be characterized. $`P(\mathrm{\Omega })`$ can be set to zero for $`\mathrm{\Omega }`$ larger than $`\kappa \omega _p`$. This assumption is reasonable below the Kosterlitz-Thouless transition. In accordance with Eq. (51), we adopt the form

$$P(\mathrm{\Omega })\frac{1}{2}((\kappa \xi ^1)\omega _p)^{1/2}\mathrm{\Theta }(\mathrm{\Omega }\omega _p\xi ^1)\mathrm{\Theta }(\kappa \omega _p\mathrm{\Omega })(\mathrm{\Omega }\omega _p\xi ^1)^{1/2},$$
(52)
Using Eq. (25) and (52), calculations can be performed analytically.

For $`eV<\mathrm{}\omega _p\xi ^1/2`$, the average current vanishes. For $`\xi ^1<2eV/\mathrm{}\omega _p<\kappa `$, the average current reads :

$$I=\frac{e}{h\mathrm{}}\frac{16\gamma ^4}{(1+\gamma ^4)^2}\frac{\kappa \omega _p}{2}(1(\kappa \xi )^1)\left(\frac{2\frac{eV}{\mathrm{}\omega _p}\xi ^1}{\kappa \xi ^1}\right)^{3/2}f(𝒜),$$
(53)
where

$`f(x)`$ $`=`$ $`{\displaystyle \frac{2}{x}}\left\{{\displaystyle \frac{1}{2}}\sqrt{1+{\displaystyle \frac{2}{x}}}\mathrm{ln}\left|{\displaystyle \frac{1+\sqrt{1+\frac{2}{x}}}{1\sqrt{1+\frac{2}{x}}}}\right|\sqrt{{\displaystyle \frac{2}{x}}1}\mathrm{arctan}\left({\displaystyle \frac{1}{\sqrt{\frac{2}{x}1}}}\right)\right\},`$ (55)
$`𝒜`$ $`=`$ $`\left({\displaystyle \frac{2eV}{\mathrm{}\omega _p}}\xi ^1\right)\left({\displaystyle \frac{\mathrm{}\omega _p}{\mathrm{\Delta }}}\right)\sqrt{1{\displaystyle \frac{4\gamma ^4}{(1+\gamma ^4)^2}}},`$ (56)
where $`𝒜`$ is a scaling variable which varies between $`0`$ (perfect contact $`\gamma =1`$) and $`1`$ ($`\gamma =0`$ and a plasma frequency $`\omega _p=|\mathrm{\Delta }|/\mathrm{}`$). Note that the function $`f`$ does not vary substantially between these two transparencies: $`f(𝒜=0)=2/3`$ and $`f(𝒜=1)0.71`$. The dependence of the current on the parameters $`\gamma `$, $`eV`$, $`\omega _p`$ and $`\xi `$ is given by the prefactor in Eq. (53).

For $`\kappa \mathrm{}\omega _p/2<eV<\mathrm{\Delta }`$,

$$I=\frac{e}{h\mathrm{}}\frac{16\gamma ^4}{(1+\gamma ^4)^2}\frac{\omega _p}{2}(\kappa \xi ^1)g(𝒜,\left(\frac{2eV}{\mathrm{}\omega _p}\xi ^1\right)(\kappa \xi ^1)^1),$$
(57)
where

$`g(x,y)`$ $`=`$ $`\left({\displaystyle \frac{y}{x}}\right)\left(\mathrm{ln}\right|1+{\displaystyle \frac{x}{2}}(1{\displaystyle \frac{1}{y}})|\mathrm{ln}|1{\displaystyle \frac{x}{2}}(1{\displaystyle \frac{1}{y}})|`$ (59)
$`\sqrt{{\displaystyle \frac{y(2+x)}{x}}}\{\mathrm{ln}|1\sqrt{{\displaystyle \frac{x}{(2+x)y}}}|\mathrm{ln}|1+\sqrt{{\displaystyle \frac{x}{(2+x)y}}}|\}2\sqrt{{\displaystyle \frac{y(2x)}{x}}}\mathrm{arctan}\left(\sqrt{{\displaystyle \frac{x}{y(2x)}}}\right)),`$
For $`\gamma =1`$, the behavior of $`I`$ is essentially the same as that of linear theory , for arbitrary $`\omega _p`$. These results are plotted in Fig. 6 (with the choice $`\omega _p=\mathrm{\Delta }/4\mathrm{}`$).

For $`eV`$ larger than $`\kappa \mathrm{}\omega _p/2`$, the curve $`I`$ versus $`eV`$ is a straight line. For less transparent barriers ($`\gamma 1`$), two behaviors have to be distinguished. First, if the plasma frequency is smaller than the gap parameter, then the scaling variable $`𝒜`$ in Eq. (56) is small and linear response is again retrieved. Second, if $`\mathrm{}\omega _p`$ is comparable to $`\mathrm{\Delta }`$, $`I`$ increases typically as $`(2eV\kappa \mathrm{}\omega _p\xi ^1/2)^{3/2}`$ for $`\xi ^1/2<2eV/\kappa \mathrm{}\omega _p<1`$. Then, for $`eV`$ larger than $`\kappa \mathrm{}\omega _p/2`$, the response is non-linear. In particular, when the transparency is low, $`I`$ increases faster than $`(2eV/\kappa \mathrm{}\omega _p\xi ^1/2)^{3/2}`$ for $`eV>\kappa \mathrm{}\omega _p/2`$. Overall, no sensible deviations from linear response are expected in the disordered limit.

## VI Conclusion

A formulation of quantum transport in normal metal superconductor junctions is developed in the case where the order parameter undergoes phase fluctuations, both in the tunnel limit and in the transparent regime. The Keldysh Green’s function was used, allowing to go beyond linear perturbation theory, provided that the subgap regime is specified (this latter assumption enabled us to neglect some contributions in the Dyson expansion). While we applied our results to the 1D Josephson junction array model of Ref. , the result is quite general as any model for the phase correlator can be employed to calculate the current and the finite frequency noise. An interesting issue is that the zero frequency noise satisfies the Schottky formula (as should be), with a current describing inelastic processes. A direct correspondence between the phase correlator and the second derivative of the noise spectral density with respect to frequency has been established. For a highly transparent interface, phase fluctuations are shown to affect only the differential conductance at low voltages: at high voltage the $`I(V)`$ characteristic is essentially that of a rigid superconductor, which can be described with the BTK model. On the other hand, for a less transparent barrier, deviations with respect to the BTK are shown to occur both at high and low voltages. This effect becomes more dramatic when the barrier is opaque.

The present system bears similarities with other situations where Andreev reflection is accompanied by the emission of a bosonic excitation. As mentioned in the introduction, Ref.  dealt with a ferromagnet whose magnons can flip the spin of one of two same spin electrons in order to contribute to the subgap conductance in a s-wave superconductor. Therefore, electrons with opposite spin contribute to the current as usual, but a new scattering channel – due to this one-magnon process – is opened, and enhances the Andreev current overall. The contrast with our calculation is that here, the current is always reduced by the excitation of phase modes. For a specified voltage bias, the more phase fluctuations, the more energetic are the modes which propagate in the Josephson junction array. On the side of the normal metal, this requires that the sum of the energies of the electron and hole has to be large. Yet the Fermi level restricts the impinging electron energy. This reduction of the available phase space for inelastic assisted Andreev reflection thus justifies a weaker current when the fluctuations are large. Also, note that the phase correlator which was used in this calculation typically appears in situations where a quantum mechanical degree of freedom is coupled to an electromagnetic environement . Overall, if one includes multi-magnon processes to the situation of Ref. , both ferromagnetic/superconductor and normal-metal/fluctuating superconductor junctions constitute illustrations of quantum dissipative systems , where the quantum process to be probed – Andreev reflection – is coupled to a physical environment (a bath of magnons/“phase oscillators”).

Potential applications to physical systems are now addressed. While typical BCS superconductors have essentially amplitude fluctuations (except close to the superconducting transition), the high-$`T_c`$ cuprates have a much lower phase stiffness and hence, are expected to show strong fluctuations of the phase of the order parameter . Here, an s-wave order parameter was assumed for simplicity, whereas these materials are mainly d-wave, but our calculation could in principle be generalized to this latter situation. The two main differences between s-wave and d-wave materials are the following: first, in a d-wave material, there are nodes in the order parameter in specific directions, which results in a finite density of states for quasiparticles near the superconductor chemical potential. Second, low energy Andreev bound states may exist at this energy. Depending on the orientation of momentum, an electron-like quasiparticle and a hole-like quasiparticle in the superconductor see an order parameter with the opposite sign . This generates a zero-bias conductance peak (ZBCP) which is a hallmark of a d-wave order parameter.

At first sight, and for ideal interfaces, the connection of between our s-wave model and a d-wave compound could be made by first choosing the antinodal direction to be orthogonal to the plane of the interface , second by restricting electron transfer along the normal of this interface. Nevertheless, the confrontation of our theoretical results with transport experiments in high-$`T_c`$ compounds may prove too difficult, as surface defects are likely to be present in such compounds. Such defects are responsible for a reduction of the order parameter in the vicinity of the junction. Moreover, the presence of Andreev bound states – also induced by the impurities at the interface – renders the present approach inadequate. Furthermore, surface roughness provokes a randomization of the tunneling directions so that the above scenario for antinodal tunneling loses its practicality.

Previous theoretical works addressing the role of phase fluctuations  on Andreev reflexion in the pseudogap phase  focussed on the mildly underdoped case: the hole density is chosen on the right hand side of the zero temperature insulator-superconductor quantum critical point (QCP). In this regime quantum fluctuations may renormalized the superfluid density but the physics turns out to be essentially that of the two-dimensional classical XY model (renormalized classical regime in the spirit of the $`n=2`$, non-linear $`\sigma `$-model). Transverse superfluid velocity fluctuations due to the motion of vortices are expected to dominate. Because the time scales for the vortex motion, as measured by THz spectroscopy  are much larger than the Andreev time $`\tau _A=\mathrm{}/\mathrm{\Delta }`$, the correlation of the phase was treated as static (on the scale of $`\tau _A`$). The phase correlator decays in space over a lengthscale $`\xi _{cl}`$ which is the classical correlation length . If tunneling occurs on an area which is much smaller than $`\xi _{cl}`$, then an Andreev signal is expected.

Our calculation (albeit in 1 dimension), focuses on quantum fluctuations. Fluctuations are included in a non-perturbative and systematic manner, regardless of the phase propagator which is used in the end. As a general trend, quantum fluctuations were found to have a tendency to reduce the Andreev reflexion signal. Note that other models of pseudogap also rely partly on quantum fluctuations , while others incorporate aspects of 1D physics (stripes) . Extensions of our calculation beyond 1D – together with surface roughness – or with competing orders could be envisioned to achieve a more realistic description of STM experiments on high-$`T_c`$ compounds.

###### Acknowledgements.

We thank G. Falci for clarifications about the phase correlator in the JJA model. Discussions with C. Bruder, D. Feinberg and F. Hekking are gratefully acknowledged.


text:
# Double Starbursts Triggered by Mergers in Hierarchical Clustering Scenarios

## 1 Introduction

The physical mechanisms that trigger and regulate star formation in galaxies are thought to be of different nature and act at different scales (e.g., Kennicutt 1998). Among them, mergers and interactions are known to play a crucial role as it has been demonstrated from observations of high emission levels of $`H_\alpha `$ (e.g., Kennicutt et al. 1987), far infrared (e.g., Solomon & Sage 1988; Sanders & Mirabel 1996) and radio continuum (e.g., Condon et al. 1982) in interacting systems. Recent observations of nearby galaxy pairs by Barton, Geller & Kenyon (2000; see also Donzelli & Pastoriza 1997) found that interactions during close encounters can be statistically correlated with enhancements of the star formation activity lasting more than $`10^8`$ yr. It has also been observed that moderately luminous starbursts in the nearby Universe may occur in disc galaxies that show only mild external perturbations originated in moderate interactions (Gallagher et al. 2001). Regarding our Galaxy, a new study of the chemical composition of the solar neighbourhood by Rocha-Pinto, Maciel & Flynn (2000) supports the possible existence of three main starburst episodes, lasting each one few Gyr, which can be associated with nearby passage of the Milky Way galaxy satellites. Remarkably, two of such starbursts show evidences of being associated to only one of those close encounters. Recent observations of high redshift ($`z`$) objects show an increase of morphologically disturbed and interacting objects with $`z`$, accompanied by an increase of the star formation activity (e.g., Bouwer, Broadhurst & Silk 1998; Driver et al. 1998; Brinchmann et al. 1998; Sawicki & Yee 1998; Le Fèvre et al. 2000). Although it is not yet clear at what extent these data are fully consistent with predictions of models based on hierarchical clustering, all of them agree in supporting a scenario in which interactions/mergers and star formation increase with $`z`$ (Flores et al. 1999; Steidel et al. 1999).

Understanding the nature of galaxies as a function of $`z`$ and the interface between galaxy formation and structure formation, requires to be able to describe how star formation proceeds in galactic objects at different stages of evolution, and how this activity is affected by the formation of the structure at larger scales. This is a complex task taking into account that as a galaxy evolves, it may react to the same physical mechanism (i.e., such as mergers/interactions) in a different way, depending on its astrophysical and dynamical properties (i.e., formation of a bulge, baryonic/dark matter dominated potential well at the centre, availability of gas, etc), and, consequently, on its previous evolutionary history.

In a hierarchical clustering scenario, galactic objects form by the aggregation of substructure. Assuming this model as the more feasible one involves the understanding of the role played by mergers and interactions as an object evolves. In particular, pre-prepared simulations of merging galaxies show that, during violent events, strong gas inflows can be triggered, feeding star formation activity at the centre (e.g., Noguchi & Ishibashi 1986; Noguchi 1991; Mihos, Bothum & Richstone 1993; Mihos 1992; Barnes & Hernquist 1991, 1996). Furthermore, Mihos & Hernquist (1994, 1996 hereafter MH96) demonstrated that during merger events, a bulge-less disc galaxy is more susceptible to growing bar instabilities that might trigger gas inflows towards its centre, even before the actual collision of the baryonic clumps occurs, than a bulge-disc system. In the first case, a double star formation burst was identified, with a first component triggered by the strong tidal fields developed during the orbital decay phase (hereafter ODP), and a second one triggered when the baryonic cores collided. The relative strength of the double bursts was shown to depend on the compactness of the bulge and on the availability of gas at the time of the merger. Conversely, a disc galaxy with a stellar bulge was found to be more stable, and only when the actual collision of the baryonic cores occurred, a starburst was produced. These authors directly link the response of the systems during a merger to the characteristics of their potential wells.

However, in these pre-prepared mergers, the great variety of phenomena involved in the complex processes of galaxy-like objects (GLOs) building-up from the field of primordial fluctuations at high $`z`$ up to $`z=0`$, cannot be followed. Such phenomena include collapse events, mergers and interactions with close neighbours, accretion of smaller objects, gas infall at scales of some hundreds of kpcs (from the environment to the halo or from the halo to the disc), gas inflow at scales of tens of kpcs (within discs and from discs to the bulge-like central concentrations), etc. All of them play a role in the assembly of galaxy-like objects, and, moreover, they interact among themselves in a non-trivial way. Pre-prepared mergers are very detailed experiments, but which only describe a particular phase of galaxy formation from initial conditions that are generally set by hand. Moreover, mergers in these experiments are taken in isolation, and not in connection with the other processes involved in galaxy formation. On the contrary, with a self-consistent approach (i.e., following the assembly of galaxy-like objects from the field of primordial fluctuations at high $`z`$ up to $`z=0`$, in the framework of a cosmological model), mergers can be studied in connection with other mechanisms, that play a fundamental role in the merger process itself, such as gas infall at some hundreds of kpc scales and gas inflows in discs. Also, one can analyze how the initial conditions that are usually considered in pre-prepared simulations (at scales of about their virial radii, i.e., some hundreds of kpcs), appear, along the evolution of the objects, in connection with the evolution of their environment at larger scales. At these larger scales, evolution is determined by the global evolution of the cosmological model. Of particular interest in this regard is to learn how bulges are assembled, and how the growth of a given bulge at a given $`z`$ might prevent it from growing at lower $`z`$s.

By using cosmological hydrodynamical simulations, Domínguez-Tenreiro, Tissera & Sáiz (1998, hereafter DTTS98) and Tissera, Sáiz & Domínguez-Tenreiro (in preparation) analyzed how gas inflows can be generated during the evolution of a typical galaxy-like object in a hierarchical clustering model, as it forms by the merger of substructure. These authors showed that the gas with non-zero angular momentum always tends to form a disc, but the survival probability of these discs depend on their stability properties. Their findings confirmed that gas inflows can be easily triggered during violent events if the axi-symmetrical character of the potential well is not assured by a central mass concentration, as for example, a dark matter halo dynamically dominating at the centre, or by the presence of a compact stellar bulge (e.g.,Toomre & Toomre 1972; Martinet 1995 and references therein). These gas inflows can be mild or catastrophic, depending on the internal properties of the disc-like systems. Catastrophic inflows can destroy the discs and are easily induced in the absence of a compact central mass concentration (DTTS98; Steinmetz & Navarro 1999). However, if stellar bulges are allowed to form, discs with observational counterparts can be formed in hierarchical clustering scenarios (Sáiz et al. 2001, hereafter Sal01).

Tissera (2000a, hereafter Tis00) studied the star formation history of galactic-like objects in relation to their merger histories within a cosmological context. This author found that minor and major mergers may trigger starbursts as powerful, regardless of their relative masses, if there is enough available gas in the system. Only the starbursts efficiencies, defined as the fraction of the available gaseous mass in the system that is actually transformed into stars, were found to depend on the relative masses of the colliding objects. The stellar masses formed and starbursts (SBs) durations measured in the simulations are in agreement with observational results (e.g., Barton et al. 2000; Glazebrook et al. 1998). Even more, Tis00 estimated a change of up to 0.6 magnitudes in the colours of the simulated galaxies due to mergers, consistent with observational results reported by Le Fèvre et al. (2000). It was also pointed out by the author that during certain merger events, double SBs were present. If these phenomena were a common physical process, then, understanding how they could be induced in a hierarchical clustering scenario is very relevant to the study of galaxy formation. The physical interpretation of several observations of galaxy properties, such as luminosity function or spatial distributions, may be affected by the possibility of inducing multiple starbursts during a single merger event (Tissera 2000b). Moreover, the formation of the stellar bulges of spiral galaxies can be directly affected by the possibility of combining two feasible formation mechanisms, such as secular evolution and collisions, during a single merger event. Observationally, there are evidences that suggest that these different processes, together with collapse at high $`z`$, could work with different relative importance, depending on the Hubble type, to manufacture stellar bulges (Ellis, Abraham & Dickinson 2001).

In this paper, we focus on the study of double starbursts within the context of galaxy formation in a cosmological framework. Contrary to pre-prepared mergers, we use fully self-consistent cosmological hydrodynamical simulations, where the distribution of merger parameters, such as orbital characteristics, the orbital energy and angular momentum, the masses of the virial haloes and baryonic clumps involved in the merger event, and the spin, internal structure, and relative orientation of the baryonic clumps that are about to merge, among others, arise naturally at each epoch as a consequence of the initial spectrum of the density fluctuation field, its normalization, and the cosmological models and its parameters. The effects of mergers and interactions on star formation can be studied in a consistent scenario, although at the expense of losing numerical resolution. In Section 2 we present the main aspects of the numerical models and in Section 3, 4 and 5, the analysis of starbursts. Section 6 summarizes the results.

## 2 Numerical Models

The simulations analyzed take into account the gravitational and hydrodynamical evolution of the matter including an algorithm to transform the cold and dense gas into stars. We use the Smooth Particle Hydrodynamical (SPH) code developed by Tissera, Lambas & Abadi (1997) coupled to the AP3M (Thomas & Couchman 1992). SPH based codes have proved to be powerful tools to study galaxy formation (e.g., Navarro & White 1994; Steinmetz & Navarro 1999; DTTS98; Tis00; Sal01).

We run four simulations, S.1, S.2, S.3 and S.4, consistent with Cold Dark Matter universes with $`\mathrm{\Omega }=1`$, $`\mathrm{\Lambda }=0`$. The normalization of the power spectrum was taken to be $`\sigma _8=0.40`$ for S.1 and S.2 and $`\sigma _8=0.67`$ for S.3 and S.4 (see Table 1). The simulated boxes have $`5h^1`$ Mpc length ($`h=0.5`$) with $`N=64^3`$ total particles. Baryonic particles represent 10 $`\%`$ of the total mass. Note that dark matter and baryonic particles have the same mass, $`M_{\mathrm{part}}=2.6\times 10^8\mathrm{M}_{}`$. The gravitational softening used is 3 kpc, and the smaller smoothing length allowed is 1.5 kpc. Simulations S.1 and S.2 share the same initial conditions (A) while S.3 and S.4 are different realizations of the power spectrum (B, C). Run S.1 corresponds to experiments S.1 in DTTS98, Sal01 and Tis00. Simulation S.2 has been run with a higher star formation efficiency in order to assess the dependence of the results on this parameter.

The star formation (SF) algorithm transforms dense and cold gas in a convergent flow into stars according to the Schmidt law (Navarro & White 1994; Tis00). We assume a certain star formation efficiency, $`c`$ (see Table 1), and a characteristic time-scale for the star formation process. The SF parameters used in S.1, S.3 and S.4 allow the transformation of gas into stars only in the very dense regions. In practice, this condition implies that gas particles have to fall in all the way down to the very central regions of GLOs before they are converted into stars (see Fig. 2). This assures the formation of a compact central stellar mass concentration, that is, a bulge, and that the transformation of gas into stars on the discs is suppressed. Conversely, the SF process in S.2 is more efficient, which implies that gas particles are transformed into stars sooner than those in S.1. This difference translates in different physical properties for the final objects, as Tissera & Domínguez-Tenreiro (1998) have shown. These authors study the mass distribution of galactic objects in S.1 and S.2 finding that, because of the higher star formation efficiency, those in S.2 tend to be dominated by stellar spheroidal components with potential wells not as deep as those of their counterparts in S.1. The results of S.2 have been included in this work to show that the conclusions we reach here are not mainly determined by the star formation parameters, implying that the mechanism at work is a general physical process.

We want to point out that, because the SF process is based on the Schmidt law, any sudden increase of the gas density directly implies a corresponding increase in the star formation activity. According to Tis00 (see also MH96), it is possible to have a large percentage of gas that satisfies the SF criterion, but that is transformed into stars at a quiescent rate. And, on the other hand, a smaller fraction of gas can be suddenly compressed, triggering a starburst.

Supernova feedback effects have not been included in these simulations. Energy feedback processes are believed to play a key role in helping to set a self-regulated star formation regimen (Silk 2001). However, these processes are not well understood and its modelling in hydrodynamical simulations is still quite controversial (Katz 1992; Navarro & White 1994; Metzler & Evrard 1994; Yepes et al. 1997; Mosconi et al. 2001). And so, as a first step, it is wise to analyze star formation separately from supernova effects. Furthermore, as pointed out by DTTS98 and Sal01, the restrictive SF model used in this work could have allowed us, in a way, to mimic the effects of supernova energy feedback. In fact, the SF algorithm we have used allows the formation of stellar bulges without quickly depleting the gas reservoirs of the objects. Sal01 have shown that the formation of compact stellar bulges provides stability to the disc-like objects, so that the disc component is able to conserve an important fraction of its angular momentum during violent events, producing systems that resemble current normal spirals. Particularly, these authors have analyzed the dynamical and structural properties of GLOs in S.1 at $`z=0`$, finding that they have observational counterparts. This is also valid for those GLOs in S.3 and S.4 that resemble disc-like objects. Hence, we are working with a set of objects that can be fairly compared to spiral galaxies at $`z=0`$.

## 3 Star Formation Histories and Mergers

Galaxy-like objects are identified at their virial radius, $`r_{200}`$ (i.e., the radius where $`\delta \rho /\rho 200`$; White & Frenk 1991) at $`z=0`$. We reject those GLOs with a massive companion within two virial radii at $`z=0`$ to avoid complications due to tidal fields originated by the presence of companions and/or the underlying overdensity, focusing on the effects produced by the assembly of each individual object through hierarchical growth. We are only going to study those GLOs with more than 250 baryonic particles inside the main object, and, generally, from $`z=1`$ to diminish as much as possible numerical resolution effects, which could be present at higher $`z`$. No other constraints have been imposed to select the analyzed GLOs.

A typical GLO is composed by a dark matter halo, a main baryonic clump and a series of small baryonic satellites. Table 2 gives the total number of dark matter ($`N_{\mathrm{dark}}`$) and baryonic ($`N_{\mathrm{bar}}`$) particles within the virial radius for each GLO at $`z=0`$. Note that dark matter haloes are very well resolved by several thousand particles. According to results from Steinmetz & White (1997), if the potential well of a galactic halo is well-reproduced, the dissipative component is forced to set onto the correct density profile, even if it is resolved by few hundred particles. If the gas density profile is adequately described, then the SF process within the GLO can be reliably followed (see also DTTS98).

The merger trees of the GLOs are reconstructed by tracking back all particles that belong to a GLO at $`z=0`$. The progenitor object is chosen as the more massive clump within this merger tree. A merger event will be defined as the whole process from the instant when the two main baryonic objects are first identified to share the same dark matter halo (this instant is the beginning of the ODP, whose redshift will hereafter be termed $`z_\mathrm{o}`$), until the instant they actually collide at a redshift $`z_\mathrm{c}`$ (once the ODP is over). Note that, in this paper, a collision means that the baryonic systems cannot be longer separated into two distinctive objects, but it does not imply any particular characteristics for the remnant. During the ODP, the minor colliding baryonic clump will be referred to as the satellite.

The corresponding star formation rate (SFR) histories of the progenitors can be constructed as the ratio between the stellar mass formed during a time-step of integration and its duration ($`\mathrm{\Delta }t=1.3\times 10^7`$ yrs for S.1 and S.2 and $`\mathrm{\Delta }t=1.2\times 10^7`$ yrs for S.3 and S.4). These SFR histories can be described as a result of two contributions, a quiescent SF and a series of starbursts (see Tis00 for details). The quiescent SF gives rise to a constant SFR or threshold; starbursts are defined to be those points of the SFR history of a given object above its threshold (see Fig. 1).

Mergers of substructure are responsible of the triggering of some of these SBs. ¿From all SFR peaks identified, in this paper we are only concerned with those that occur within merger events. We have studied 29 merger events related to 43 starbursts of which 14 are triggered during orbital decay phases. Once the starbursts are isolated, it is direct to define their maximum SFR, $`\sigma _{\mathrm{star}}`$, the total stellar mass formed, $`M_{\mathrm{burst}}`$, and their duration, $`\tau _{\mathrm{burst}}`$. These quantities determine the characteristics of each starburst, together with the gas richness of the system, $`M_{\mathrm{star}}/M_{\mathrm{bar}}`$, and the virial mass ratio of the colliding objects, $`M_{\mathrm{sat}}/M_{\mathrm{pro}}`$, at the time of the merger. $`M_{\mathrm{star}}`$ and $`M_{\mathrm{bar}}`$ are the stellar and baryonic (i.e., gas plus stars) content of the merging systems within $`r_{200}`$; $`M_{\mathrm{sat}}`$ and $`M_{\mathrm{pro}}`$ are the virial masses of the satellite and progenitor, respectively. Table 2 gives their values for the mergers analyzed in this work.

## 4 Double Starbursts

During some merger events, two starbursts can be identified in the SFR histories of the GLOs. We classify as a double SB those two well-defined SFR peaks that could be directly related to the same merger process. In other words, those two bursts that appear in the time interval from $`z_\mathrm{o}`$ (i.e., the instant when the satellite enters the virial radius of the progenitor) to the actual collision of the baryonic cores at $`z_\mathrm{c}`$. Only those bursts that strictly satisfy this condition were taken as probably induced by the same merger. Whenever a doubt exists, we choose to eliminate the merger from our statistics so that it does not contribute as either, a single burst or a double one. This situation, actually, occurs once. As an example, we show in Fig. 1 the SFR for typical GLOs that experience (a) a double SB (GLO 3 in S.1, Fig. 1a) and (b) a single one (GLO 3 in S.3, second merger, Fig. 1b) as a function of lookback time ($`\tau (z)=1(1+z)^{3/2}`$). We also plot in each of these Figures the distances from the progenitor to the satellite in the time interval from $`z_\mathrm{o}`$ to $`z_\mathrm{c}`$<sup>1</sup><sup>1</sup>1 It is worth mentioning that in some merger events, during the orbiting of the satellite around the progenitor, more than one SBs can be detected. These cases are not studied in this paper.. In Fig. 1a, we see that the first burst occurs during the orbital decay phase, while the second one is triggered during the collision of the baryonic cores. In Fig. 1b, only at the impact of the two baryonic clumps an enhancement of the SF activity is detected. Note that when the satellite enters the virial radius of the progenitor, in both cases (a, b), the progenitors are equally gas-rich as can be appreciated from Table 2.

For each double SB, we measured the $`\tau _{\mathrm{burst}}`$, $`M_{\mathrm{burst}}`$, $`\sigma _{\mathrm{star}}`$ and the time interval $`\mathrm{\Delta }t_{\mathrm{burst}}`$ between the primary ($`\sigma _{\mathrm{star}}^1`$) and the secondary ($`\sigma _{\mathrm{star}}^2`$) burst maxima, ($`z_{\mathrm{burst},1}>z_{\mathrm{burst},2}`$; $`\mathrm{\Delta }t_{\mathrm{burst}}`$ was measured from the last SFR value of the primary burst to the first point of the secondary one, above the thresholds defining the bursts). The temporal separations between the bursts, $`\mathrm{\Delta }t_{\mathrm{burst}}5\times 10^85\times 10^9\mathrm{yrs}`$, are consistent with the collapse times of satellites estimated from dynamical friction considerations (Navarro, Frenk & White 1995; Hernquist 1989). Table 2 summarizes the results (secondary components of double bursts are indicated by the letter D). Last column on Table 2 gives the redshifts ($`z`$) at which each starburst starts.

It is clear that if there is enough available gas in the system, a starburst will be triggered when the progenitor and the satellite collide at $`z_\mathrm{c}`$. Regarding double SB events, our concern now is to prove that the SBs detected in the progenitors during some ODPs are produced by gas inflows driven as the satellite gets closer, as predicted by pre-prepared simulations (e.g., Hernquist 1989; MH96). Also, we want now to assess the differences, prior to the merger event, of the internal structure of GLOs involved in double and single SBs, as well as the evolution of these differences, in relation with their SFR histories, during and immediately after the merger event. To this end, we have studied and compared the properties of our GLO sample at the merger events reported in Table 2. For each of these mergers, the variations of the baryonic mass and angular momentum distribution at the central regions of their respective progenitors have been followed and quantified, in the time interval between the beginning of their ODPs at $`z_\mathrm{o}`$ and a reference time at $`z_\mathrm{f}`$<sup>2</sup><sup>2</sup>2 For mergers producing double starbursts, $`z_\mathrm{f}`$ has been taken to be the moment when the ODP ends, $`z_{\mathrm{o},\mathrm{f}}`$, because $`[z_\mathrm{o},z_{\mathrm{o},\mathrm{f}}]`$ is the time interval within which the primary component of double starbursts is detected. For mergers that do not produce double starbursts, we have taken $`z_\mathrm{f}=z_\mathrm{c}`$, that is, the time when the satellite and the progenitor actually collide., using different estimators.

As a first step, the progenitor particle number (actually, the baryon mass) inside a given radius, $`N(<r)`$, has been measured at different redshifts in the $`[z_\mathrm{o},z_\mathrm{f}]`$ time interval. The general result is that gas inflows are detected during this time interval at the central regions of progenitors that experience double starbursts. In contrast with this result, no gas inflows are detected at any stage at the central regions of progenitors involved in merger events with single starbursts. To illustrate these results, we have chosen as typical examples the mergers occurring in GLO 2 from S.4 (a double SB event at $`z1.20`$) and in GLO 5 from S.1 (the latest one at $`z0.19`$, a single SB event; see Table 2), but the behaviour is similar for the other merger events. In Fig. 2a we plot $`N(<r)`$ for the primary component of the double SB in GLO 2 from $`z_\mathrm{o}=1.81`$ to $`z_\mathrm{f}=1.28`$. Note how gas particles are displaced from the outer to the inner radius as the time goes on. Part of the particles that get into the 10 kpc $`<r<`$ 30 kpc shell are transported all the way down to the very centre ($`r<`$ 4 kpc). This increase in the gas density immediately translates into an increase of the star formation activity, as previously discussed. The inward displacement of these particles is due to a loss of their spin angular momentum content as the satellite gets closer (see DTTS98). GLO 5 has assembled most of its mass at higher $`z`$, so that when the last merger occurs, its gaseous disc and stellar bulge are clearly at place. Fig. 2b shows $`N(<r)`$ between $`z_\mathrm{o}`$ and $`z_\mathrm{f}`$ for the last merger. As can be seen, the changes in the baryonic mass distribution are much less important in this case than for GLO 2. There is almost no increase in the baryonic mass within 15 kpc or the inner 4 kpc, implying that no significant gas inflows have been triggered in the central region during the orbital decay phase of this satellite.

In the other six plots of Fig. 2, we show the angular momentum component per unit mass, $`j_{z,i}`$, along the direction of the total angular momentum of the main baryonic systems, for each gas and star particle in GLO 2 (left panels) and GLO 5 (right panels), versus particle radial distance, $`r_i`$, at three stages of evolution within the interval of interest $`[z_\mathrm{o},z_\mathrm{f}]`$. These plots give a clear picture of the internal structure of both GLOs during such interval. These Figures show that when GLOs have well-formed stellar bulges, the disc components are not significantly affected during the ODP. Conversely, those GLOs which experience early gas inflows do so at the expense of the gas component on the discs. These gas inflows feed the stellar bulges. Fig. 2 also suggests that in the latter cases, GLOs are in early stages of evolution.

To go further into the quantification of the differences between the properties of the GLOs at a double or single SB event, we use results from bulge-disc decomposition of their projected mass surface density (Scannapieco et al., in preparation). Following Sal01, the projected mass surface density has been decomposed by using a Sérsic law for the bulge component and an exponential law for the disc component. Sal01 have shown that, at $`z=0`$, this procedure yields scalelengths and dynamical parameters which are in very good agreement with recent observations of spiral galaxies (e.g., Courteau 1996). Here, the aim of this decomposition is to provide structural parameters which allow a convenient quantification of the mass distributions of progenitors at $`z_\mathrm{o}`$ and of their changes between $`z_\mathrm{o}`$ and $`z_\mathrm{f}`$. The analysis has been performed for GLOs in S.1 and S.4, for which a suitable number of outputs to follow their evolution in time are available. The general result is that the mean values of the masses of the bulges at $`z_\mathrm{o}`$, as well as their variations between $`z_\mathrm{o}`$ and $`z_{\mathrm{o},\mathrm{f}}`$, show clear differences for GLOs involved in double or single starburst events. At $`z_\mathrm{o}`$, the former show smaller bulges than the later on average (i.e., $`50\%`$ smaller), but, after the starbursts that take place within the $`[z_\mathrm{o},z_{\mathrm{o},\mathrm{f}}]`$ interval, their bulges have grown by $`47\pm 19\%`$ on average. Meanwhile, those GLOs with no important SF during an equivalent time-interval, show a change in the bulge mass of $`5\pm 2\%`$<sup>3</sup><sup>3</sup>3 The errors correspond to standard deviations calculated by using the re-sampling bootstrap technique with 500 random samples.. These results confirm that systems experiencing early gas inflows have smaller bulges which are fed by gas inflows during the ODPs.

## 5 Starburst Properties

Taking into account previous works on disc formation, there are evidences that support the idea that the stability of the disc structures plays a relevant role in the triggering of gas inflows when the systems are subject to tidal fields. This stability could be assured by a compact stellar bulge when the dark matter (DM) halo is not dynamically dominating at the centre (MH96; Barnes & Hernquist 1996). Then, a first step towards the understanding of the physical origin of the induced SBs in our simulations, is to analyze any possible relation between the formation of the stellar bulge in the progenitor and the burst properties.

To this end, our first task is to define the epoch of bulge formation. This is a controversial point since bulges can continuously grow by mergers, bar-driven infall or other mechanisms (Combes 2000). The simplest possibility to define the presence of a bulge is to adopt a threshold criterion relative to the bulge mass. However, we want to employ a criterion which allows the comparison among bulges of different virial masses. Consequently, due to the correlation among bulge masses and virial masses, a threshold criterion relative to normalized bulge masses should be used instead. Moreover, the normalizing factor has also to be sensitive to the history of evolution of GLOs (i.e., their respective merger trees and interaction events), because this history also affects the SF history of each GLO. We found that the total stellar mass in the GLO at $`z=0`$ is sensitive to both the virial mass of the object and its history of evolution. Then, we define the ratio between the stellar mass formed in the progenitor at a given $`z`$, $`M_{\mathrm{star}}^z`$, and the total final stellar mass, $`M_{\mathrm{star}}^{z=0}`$, of the GLO, as an estimator of the stellar core (i.e., bulge) formation (see Fig. 2 for a general picture of baryon distributions). The monitoring of this ratio along the history of evolution of a GLO allows us to adopt a clean criterion for the formation of bulges.

Now, if our hypothesis about the ability of well-formed bulges to provide stability to a bulge-disc system is correct, then there should be a correlation between the $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$ ratios, measured at a redshift $`z=z_\mathrm{o}`$, when the satellite enters the virial radius of the parent object, and the properties of the starbursts. For each merger in Table 2, we have calculated $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$ and the relative strengths of the two bursts, $`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^2`$. The results are plotted in Fig. 3. As can be seen from this Figure, there is a clear segregation which indicates that when a single burst ($`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^2=0`$) is associated with a merger, GLOs tend to have larger $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$ values. Conversely, when a double burst is found, all of them have $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}0.40`$. Our results show that as the stellar bulge forms, the relative strength relation reverses ($`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^2<1`$) until the first burst disappears when the stellar bulge is fully present ($`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^2=0`$). The ratio $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}0.40`$ seems to indicate when the stellar bulges have grown large enough to prevent even mild inflows.

Hence, it seems plausible that, depending on the structural properties of the GLOs, the merger with a satellite may induce two starbursts whose properties would also depend on the internal structure of the galactic objects. Our findings also agree with Hernquist (1989) results, which suggest that tidal triggering of starbursts may be a natural way of explaining the discrepancy found between the inferred starburst time-scales ($`10^8\mathrm{yr}`$) and the dynamical time-scales of the objects ($`10^9\mathrm{yr}`$) .

### 5.1 Dependence on the potential well

In order to intend to relate the characteristics of the potential well in the very central regions of GLOs with the presence of a well-formed bulge as measured by $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$, we estimate the circular velocity at $`r=`$ 2 kpc, $`\upsilon _k`$, for the dark matter ($`k=`$ DM) and gaseous ($`k=`$ gas) components and the total mass ($`k=`$ tot) at $`z_\mathrm{o}`$. This is a simple but straightforward way of measuring the mass concentration. For this purpose, we calculate the circular velocity curves defined as $`V(r)^2=GM(r)/r`$, for the three mass components at primary components and single bursts. We use velocities instead of density profiles since the former are defined by the integrated mass which is less affected by numerical noise. Note also that for this calculation we are not including the gravitational softening since we are interested in measuring the actual mass distribution. However, the correct circular velocity curves should be estimated by taking into account such softening as shown by Tissera & Domínguez-Tenreiro (1998) and Sal01.

Figs. 4 (a,b,c) show $`\upsilon _k`$ for the total mass, DM and gaseous components, respectively, against the ratio between the heights of the starbursts $`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^2`$ (recall that for single bursts, $`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^2=0`$) for bursts in S.1 and S.4. The segregation is quite clear, indicating that only when the systems have less concentrated potential wells measured by lower $`\upsilon _k`$ values, gas inflows are induced, triggering starbursts during the orbital decay period. Note that all mass components show a similar trend. As it can be appreciated from Fig. 3, systems with $`\sigma _{\mathrm{star}}^1/\sigma _{\mathrm{star}}^20`$ have smaller stellar bulges. Because of the new gas infall, the stellar bulges grow as the gas is transformed into stars at the centre. When GLOs have well-formed bulges, they have already pulled in enough mass to assure the stability of the system which prevent further important inflows.

### 5.2 Starburst efficiencies

¿From Fig. 3 we showed that double bursts occur only when $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}0.40`$. It was found by Tis00 that the efficiency in the transformation of gas into stars during SBs correlates with $`M_{\mathrm{sat}}/M_{\mathrm{pro}}`$ in the sense that massive mergers ($`M_{\mathrm{sat}}/M_{\mathrm{pro}}0.40`$) seem to be more efficient at converting gas into stars. Here we explore the possibility that the SB efficiencies (defined as the ratio between the stellar mass formed during a burst and the total gas mass available in the system at the time of the merger, $`M_{\mathrm{burst}}/M_{\mathrm{gas}}`$) of the bursts triggered during the ODP may be different from those of SBs produced at the actual collision of the baryonic clumps. In double SBs we do not consider the secondary components, since they form from the left-over gas in the progenitor after the primary one and the gas component associated to the satellite. They are also difficult to study since we would need to temporary resolve the merger with a higher number of outputs and higher numerical resolution.

Fig. 5 shows $`M_{\mathrm{burst}}/M_{\mathrm{gas}}`$ vs $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$. Among GLOs in the same simulation, the larger values are mainly recorded for the smaller $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$. This can be clearly appreciated from Table 3 which shows the mean values of $`M_{\mathrm{burst}}/M_{\mathrm{gas}}`$ for $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$ smaller or greater than 0.40. Note that mean starburst efficiencies are higher for GLOs in S.2 simulation.

Hence, we find a clear trend for SBs induced during the ODP to be more efficient at converting gas into stars than those produced by the actual collision of the baryonic cores. This result may have strong implications for galaxy formation.

### 5.3 Dependence on redshift

In Figs. 6a and 6b we plot $`M_{\mathrm{sat}}/M_{\mathrm{pro}}`$ and $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$, respectively, as a function of the redshift when the starbursts start ($`z_{\mathrm{burst}}`$ on Table 2). We include single bursts and primary components of double SBs (features with arrows). As can be seen form Fig. 6a, there is a clear trend to have major mergers (i.e., similar mass objects) at high $`z`$ as expected in a hierarchical clustering model (e.g., Lacey & Cole 1993). We also see that double SBs tend to occur at higher $`z`$, and consequently, to be linked to more important merger events, with $`M_{\mathrm{sat}}/M_{\mathrm{pro}}0.20`$ for all double events. In Fig. 6b, $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$ as a function of $`z_{\mathrm{burst}}`$ shows how, as the progenitors build up their stellar bulges as they evolve, the probability of experiencing bursts induced via tidal forces diminishes. Note that, in these models, all double events occur at $`z>0.4`$.

¿From the combination of the information coming out from Fig. 6, we note that $`61\%`$ of bursts induced during the ODP are triggered by minor mergers ($`M_{\mathrm{sat}}/M_{\mathrm{pro}}0.40`$). Even more, if we consider only SBs in S.3 and S.4, which constitute an homogeneous sample (i.e., different random phases of the power spectrum run with the same model parameters), the trend is very similar: $`66\%`$ of the SBs detected during the ODP are related to minor mergers. This tendency suggests that the relative mass of the colliding objects is not the determining mechanism in the triggering of the early gas inflows. On the other hand, none of the systems with well-formed stellar bulges, $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}>0.40`$, have starbursts induced during the ODPs (Fig. 3).

¿From these results, we conclude that hierarchical clustering is determining the global trend for massive mergers to occur at higher $`z`$, when the systems are more likely to be unstable, and, consequently, more susceptible to experiencing gas inflows and double SB events. We predict that this trend to have SBs induced via tidal forces at higher $`z`$, combined with their higher efficiency in transforming gas into stars, lead to a global star formation efficiency that would increase with $`z`$. Observational analysis of the star formation efficiency with $`z`$ should test this prediction.

### 5.4 Bulge formation

Three different formation scenarios have been proposed for bulges of spiral galaxies: dissipative collapse (Eggen, Lynden-Bell & Sandage 1962), mergers of previously formed systems (e.g., Kauffmann, White & Guiderdoni 1993) and secular evolution (see Combes 2000 for a review). It has been also suggested that the three of them may have played a role, but with different relative efficiencies depending on the Hubble type (e.g., Combes 2000; Ellis et al. 2001). Our results support such hypotheses. Moreover, they show that a given merger event could contribute to the formation of stellar bulges through two different processes: a) gas inflows triggered during interactions at the ODPs, and b) collision of baryonic cores, once the ODPs are over. The first contribute to the secular evolution of bulges and have been found to be more important at higher $`z`$, as previously discussed (see Figs. 6; note that these inflows are episodic). Together with this secular component, there is a collisional contribution produced during the actual mergers of the main baryonic clumps, so that after a merger event, the stellar population of the GLO involved in the merger consists of the contribution of both components, plus the old stars. The relative importance of each of these stellar components would depend on the previous evolutionary paths of the objects and their physical internal characteristics at the time of each merger, particularly on those of their potential wells (see subsection 5.1).

According to this scenario, the stellar bulges will have stellar populations with different ages. For the particular models analyzed in this paper, $`50\%`$ of the stars in the bulge are formed at $`z0.5`$. If we assume that, morphologically, the bulge component is at place when $`50\%`$ of its total stellar mass is formed, then $`z0.5`$ will be its average age in these CDM scenarios. However, if the older stellar populations are taken to be age indicators, then the time of formation is displaced to much higher $`z`$.

The time of formation of the bulge component could be affected by the cosmological model as well as by the treatment of supernova effects. In particular, the lower normalization parameters ($`\sigma _8`$) adopted for these simulations contribute to produce a late formation of bulges. In a higher normalization scenario, mergers would occur at higher $`z`$, with an earlier formation of central mass concentrations. Note also that higher numerical resolution can push the age towards higher redshift. Hence, our age indicators should be taken as lower limits. Nevertheless, the physical processes and, principally, the fact that stellar bulges are found to be formed naturally as the product of collapse, collisions and secular evolution, remain valid. Although previous works proposed such mechanisms to work together in the manufacturing of stellar bulges, the results of the simulations reported here constitute the first proof that these mechanisms are actually at work in a cosmological framework and, also, they give us a preliminary information on how they do take place. The study of bulge formation by means of numerical simulations is very convenient, since it can provide hints on the physical processes involved and their interplay, while in semi-analytical models they have to be modelled following recipes known in advance (e.g., van den Bosch 2001).

Observationally, it has been recently shown by Ellis et al. (2001) that at $`z0.6`$ the bulges of spiral galaxies are bluer than ellipticals, suggesting that the stellar population has been rejuvenated, specially at intermediate $`z`$. This has been suggested to be in contradiction with predictions of some hierarchical clustering scenarios made by semi-analytical models, where spiral bulges and ellipticals are found to have very similar characteristics, as, in these models, their stellar populations have formed mostly at collapse time. However, our results show that hierarchical clustering scenarios can be reconciled with these observations, since the bulge stellar populations formed through mergers events, might be the rejuvenating agent after core collapse.

## 6 Discussion and Conclusions

The main advantage of the method reported in this paper to study mergers, over pre-prepared simulations, is that we analyze mergers that occur as the result of the self-consistent formation of the structure within a cosmological model. Hence, the physical properties of the colliding objects and the parameters of the encounters are those that correspond to galactic objects evolving in the adopted framework, without any external manipulation.

In particular, the processes of gas infall at scales of hundreds of kpcs and gas inflows in discs, at smaller scales, appear in a natural way in our simulations. Both of these processes are a consequence of evolution at larger scales, from initial conditions that sample a given spectrum of primordial fluctuations (CDM in this case) at very high $`z`$. The importance of these processes in the SFR history of GLOs lies in the role they play, not only as suppliers of gas to be turned into stars in GLOs, but also in the fact that they provide the baryon mass to form a central compact concentration, that, later on, makes it difficult any further gas inflow. It is this possibility of analyzing different phenomena involved in galaxy formation, and the dynamical and hydrodynamical connections among them, that renders attractive the self-consistent approach to study mergers reported in this article.

Our results suggest that the internal structure of the GLOs affects its SFR history, as reported by MH96 among others. Because a disc structure is a natural attractor for the dissipative component, its stability is crucial for its later evolution. The possible formation of radial inflows which can re-distribute the gas within the central regions and accumulate it in very short time-scales, directly affects the SFR in the systems.

In agreement with previous works, we find that some mergers can induce star formation in two stages, depending on the internal mass distributions of the parent galaxy. The second burst is simultaneous with the collision of the baryonic cores, while the first one occurs during the orbital decay phase of the satellite. We also found that the starburst efficiency of the primary bursts tends to be larger than that of the secondary ones, implying that induced starbursts during the orbital decay phase can be more efficient at transforming gas into stars than those produced in the actual collision of the baryonic cores. Therefore, although the rate of gas cooling and accretion onto the main baryonic clump, at scales of say $``$ 100 kpc, can be mainly determined by the dark matter halo (White & Frenk 1991), how this gas distributes itself and responses to external and/or internal factors within the central GLO regions is a small scale process and it is affected by its particular evolutionary history (DTTS98; Sal01).

A specific result of our method of analysis is the finding that a self-regulated process among gas, stars and dark matter components is settled. This self-regulation is a basic clue to understand bulge formation that can only be studied in a self-consistent scenario. Our results show that single SBs tend to occur in systems with larger $`M_{\mathrm{star}}^z/M_{\mathrm{star}}^{z=0}`$, which is pointing out the fact that these systems have suffered gas inflows in previous phases of evolution which have fed their bulges. As a consequence, the dark matter has also got more concentrated as the response to the baryon concentration (see Tissera & Domínguez-Tenreiro 1998). Conversely, systems with shallower potential wells tend to be in the process of forming the bulge. During this process, the whole potential well gets deeper. The three components, gas, stars and DM, work together. It is the response to their coupled evolution that can be correlated with the effects a merger event may cause on the dynamics of the systems and, as a consequence, on their star formation histories. Other characteristics of mergers, such as orbital orientation, may also have influence on the response of the systems as reported by Barnes & Hernquist (1996). Their detailed study within a hierarchical scenario will be carried out in a separate paper.

We found that in a hierarchical clustering scenario, there is a natural evolutionary sequence: mergers occurring at higher $`z`$ when the systems are not that concentrated and stellar bulges not well-formed, lead to gas inflows that trigger SBs and feed the first bulges. As a consequence, the total potential wells concentrate. When the objects reach lower $`z`$, they tend to be more stable because of their important stellar bulges and concentrated potential wells that prevent further gas inflows. At this stage, mergers also tend to be minor. This evolutionary sequence suggests that starbursts induced during the ODP are likely to be more common at higher redshifts; and because of their higher efficiency in transforming gas into stars, we predict that their impact on galaxy formation could be very important, particularly, on the formation of spiral bulges.

Finally, our results suggest that hierarchical clustering scenarios can be reconciled with recent observational results of high redshift spiral bulges and ellipticals in the field, since, in such scenarios, stellar bulges would form by the contribution of collapse, collisions and secular evolution, each one with a relative importance determined by the evolutionary path of each system.

We are grateful to the University of Oxford (United Kingdom), Imperial College of Science, Technology and Medicine (United Kingdom), Observatorio Astronómico de Córdoba (Argentina) and Centro de Computación Científica de la Universidad Autónoma de Madrid (Spain) for providing the computational support for this work. This work was partially supported by MCyT (Spain), through grant AYA-0973, Consejo Nacional de Ciencia y Tecnología (Argentina) and Agencia de Promoción Científica y Técnica (Argentina).

## Figure captions