Kosterlitz-Thouless Devillard

Josephson junctions chains with long-range interactions, phase slip proliferation versus

Keywords: numbers: 73.23.-b, 74.45.+c, 74.81.Fa

The role of long range badly screened Coulomb interactions in a one-dimensional chain of Josephson junctions is studied. Correlation functions for the phase correlator are obtained as a function of the Josephson coupling energy, the short range part of Coulomb repulsion and its long range component. Though quasi-long range order is no longer possible and the usual Kosterlitz-Thouless transition no longer exists, there are remnants of it. As an application, we calculate the I -V curves for Andreev reflexion when a normal metal is placed in contact with the chain. Formally, there is always an offset voltage V0 below which no current can flow, however, in some regimes V0 can be negligible. Contrary to what happens without long-range interactions, the Andreev current, as a function of applied voltage, increases faster than any power law. Signatures of long range interactions and phase slips appear in the I -V curves. Possible application for quasi one-dimensional thin superconducting wires is outlined.

I. INTRODUCTION

Josephson junction arrays (JJA) have been widely studied, both theoretically and experimentally 1 . They show a number of phenomena, vortex physics 2 , quantum phase transitions. They are also a paradigm; a large variety of physical systems can be modelled using Josephson junctions arrays, for example, one-dimensional thin superconducting wires 3,5 and granular superconductors. Nowadays, they are used to realize protected qbits 4 or for metrological applications, for establishing current standards through frequency measurements of Shapiro steps 6 .

Here, we shall restrict ourselves to the 1-dimensional case. Furthermore, we shall neglect dissipation, which is crucial in some quantum wires. No coupling to a dissipative environment is considered. For normal regular junction chains, the essential physics is the competition between Josephson coupling and local Coulomb interactions. In one dimension, true long range order can not exist, the correlations of the superconducting order parameter decay at best as a power law. When local Coulomb interactions are increased, the power law is modified, and eventually a Berezinskii-Kosterlitz-Thouless (BKT) transition occurs. Beyond the transition, when Coulomb interaction dominates, the decay of correlations becomes exponential. This is due to the unbinding of phase slips, whose existence have been shown directly in recent experiments 6 .

Topological defects have always played a central role, from BKT transitions in magnets, roughening of interfaces, disordered Luttinger liquids and some stripe models for high-T c superconductors.

We would like in this paper to study the role of long range Coulomb interactions, hereafter denoted as LRI, which have not been taken into account so far in one-dimensional JJ chains. One motivation is that, in one-dimensional clean superconducting wires, the long range part of the Coulomb interactions is poorly screened. For example, the plasma mode frequency goes to zero for small frequency 7 . This has motivated some recent studies 8 . Clearly, in these systems, fluctuations of the phase of the order parameter have to be taken into account. The Hamiltonian looks very much like the one of a one-dimensional Josephson array, with a long range Coulomb term in addition.

Long range interactions, hereafter denoted by LRI, were studied a long time ago at the spin-wave level for electronic systems 9,[START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF] . They also appear naturally in modelling underdoped high-T c superconductors, where phase fluctuations of the order parameter are important [START_REF] De Palo | [END_REF]12 . Influence of LRI have also been studied recently 13 theoretically in two and three dimensions. Our goal is to obtain the phase correlator of the superconducting order parameter as a function of time and space. This quantity governs numerous physical properties. In this paper, we shall give only one application. Namely, we calculate the average Andreev current from a normal metal to a JJ chain. Because LRI are normally screened in fabricated JJ chains, the model discussed here is not meant to be a realistic model for a JJ chain but can be relevant to thin 1-dimensional superconducting wires. However, calculations developped in this paper can be adapted to take into account the screening of the Coulomb interaction between islands and make predictions about Andreev current in real JJ chains.

After the introduction, the body of the paper is organized in four parts.

In part two, we first describe the model of JJ chain we use, and its connection to high-T c materials. It is convenient to go to imaginary time. When the Josephson coupling dominates over the Coulomb repulsion, if phase slips are neglected, the modes are spin waves. The long range interactions cause the decay of correlations to change from those of a Luttinger liquid to a more pronounced decay 14 .

In the next section, we take into account phase slips. The Hamiltonian decouples into a part describing spin-waves (long wavelength excitations) and another part describing phase slips (rapid variations of the phase). In the last sections, we discuss the interaction between phase slips. For distances much smaller than a characteristic distance, they interact logarithmically. Beyond this distance, the interaction behaves only as the square root of the logarithm of the distance 15 .

In part three, we examine how the KT real space renormalization group (RG) analysis is modified by LRI. There, as soon as LRI are present, strictly speaking, no quasi long range order is possible anymore. However, the initial long range interactions can be very weak. If we iterate the RG flow, it will take some time for LRI to alter the behavior of the usual KT flow. The LRI will show up only after a certain crossover "time", which we calculate and corresponds to some temperature T * . On the other hand, RG has to be stopped at some cutoff "time", /k B T . Therefore, in order to see the influence of LRI, the temperature must be lower than T * . The dependence of T * on the values of parameters such as the Coulomb energy, Josephson coupling and strength of LRI is explicited.

Part four is devoted to the calculation of the phase correlation function. Analytical continuation is required to obtain the correlator in real time. This part is rather technical; some aspects of the calculations are relegated to appendices. The main result is the expression of the phase correlator, which differs from the usual Ornstein-Zernike expression.

In the last part, we discuss one physical situation. The Andreev current between a metal and a JJ chain, as a function of applied voltage is calculated. We compare the curves obtained in three cases. The first case corresponds to only on-site Coulomb interaction. In the second case, a Coulomb interaction between nearest-neighbor islands is added. This model is often used in the literature. Finally, in the third case, LRI are added. For relatively low voltages, the Andreev current remains always smaller with LRI. We focuse on the case where, without LRI, we would be in the disordered phase. There is a threshold voltage V 0 below which no Andreev current can flow. This remains true in the presence of LRI. However, the current versus voltage (I -V ) curves are different from the case without LRI. Without LRI, the Andreev current starts as (V -V 0 ) 3/2 , whereas with LRI, it starts to increase faster than any power law. This can be viewed as a signature of LRI and could be tested in experiments.

In the conclusion, we discuss two possible applications, one for realistic fabricated JJ chains and the other for one-dimensional superconducting wires. For JJ chains, the interaction between island charges is screened but the methods used in this paper can be adapted to make at least qualitative predictions. For 1d wires, owing to the strong dependence of a phase slip fugacity with temperature and wire width, and also for other reasons detailed below, it is not obvious that LRI can be detected experimentally. Maybe they could for very long wires.

II. GENERALIZED VILLAIN MODEL AND EFFECTIVE ACTION WITH LONG RANGE, XY MODEL IN TWO DIMENSIONS

A. Imaginary time action

We would like to study the influence of poorly screened Coulomb interactions in 1 dimension. One possibility, which has been widely used some time ago, mainly in dimensions two and three, would be to take the Hamiltonian for a BCS superconductor [START_REF] De Palo | [END_REF]12 and to introduce Hubbard-Stratonovich transformations in the particle-hole and particle-particle channels in order to make the usual gap ∆ = |∆|e iϕ appear. The amplitude fluctuations are integrated out, so are the phase of the density fluctuations, leaving us with a phase action only. This procedure has been used in Ref. 8, except that they stayed at the RPA level and did not consider topological excitations for the phase ϕ of the order parameter. Another way would be to start with a one-dimensional array of Josephson junctions, modelled by an ideal Josephson junction (JJ) and a capacity C J in parallel. In addition, on each island, the capacity to the ground has C g to be taken into account. See Fig. 1.

The Hamiltonian is

H JJA = i J cos(θ i+1 -θ i ) + 1 2 i,j n i U i,j n j + A 2C i,j q i q j |i -j| , (1)
with q i = 2en i and the interaction potential U i,j = 2e 2 C -1 i,j and C -1 i,j is the capacitance matrix. The last term of Eq. (1) represents LRI, A is a constant and C a capacity. Usually, only the capacity of one junction C J and the capacity from an island to the ground C g are taken into account, the capacity matrix C i,j has the form

C i,j = C g + 2C J if i = j, C i,j = -C J if i = j ± 1,
and C i,j = 0 otherwise. It is reasonable to approximate the interaction potential 32 by

C g C g C g C J C J C J i Q j Q normal metal Island i Island j
Josephson junctions chain long range interaction the form

U i,j = 2 e 2 C J λ exp - |i -j| λ , (2)
with λ = C J /C g . The screening length is of the order aλ, with a the interjunction spacing. Typically, for onedimensional fabricated chains 34 C J ≫ C g but still λ is of order 10 to 30, not more. LRI thus do not appear to be completely realistic for fabricated JJ chains, but for thin superconducting wires, the Coulomb interaction is badly screened [START_REF] Giamarchi | Quantum Physics in One Dimension[END_REF]28 and can be longer range. We shall thus take a Hamiltonian of the form

H JJA = i J cos(θ i+1 -θ i ) + 1 2 i,j q 2 i 2C + A 2C i,j q i q j |i -j| , (3)
with C -1 typically an effective capacitance to the ground. We can also see this Hamiltonian as a 1d analog of the following Hamiltonian for superconductors with fluctuations of the phase of the order parameter

H S = i 2 N 0 S 2m * cos(θ i+1 -θ i ) + i 1 2 m * v 2 N 0 S ρ 2 i + i,j πe 2 ǫ ρ i ρ j |r i -r j | . (4)
In this case, N 0 S is the superfluid density in 3d. ρ is the charge density, θ will denote the phase of the order parameter. Its conjugate variable is Π, which is proportional to the charge q. When phase slips are neglected, the procedure to obtain an effective action for θ is standard and we just recall it here for completeness.

Writing q i = 2e n i with n i integers, θ i and n i are conjugate variables, [n i , θ i ] = δ i,j ,

H = i,j -J cos(θ i -θ j ) + 4E 0 n 2 i + 4AE 0 i,j n i n j |i -j| , (5)
with E 0 = e 2 2C the charging energy for nearest neighbor. Following Ref. 16, T denoting the temperature and β = 1/(k B T ), the partition function is calculated in imaginary time, by discretizing in time steps of size / √ 2E 0 J. The system without LRI would be equivalent to a usual XY model with isotropic coupling. In the following, θ i,j will no longer be operators but numbers; index i from 1 to N ≡ L/a denotes the position in real space and the second index j from 1 to N ′ ≡ β √ 2E 0 J the position on the imaginary time axis. Making = 1, the partition function reads

Z = {θ} exp - i,j J 2E 0 [cos(θ i+1,j -θ i,j ) + cos(θ i,j+1 -θ i,j )] , (6)
where the sum is on all configurations {θ}. Now, in the spin-wave approximation, the cosines are expanded and the Lagrangean L and action S take the form

L = - N i=1 1 2 J(θ i+1,τ -θ i,τ ) 2 + 4E 0 n i (τ ) 2 + 4AE 0 i i ′ n i n j |i -j| + i n i (τ )∂ τ θ i (τ ), (7)
S = -i β 0 L dτ. (8)
Applying a Fourier transform and integrating out the Π q = n q field,

S ef f = β 0 i 1 2 J θ i+1 (τ ) -θ i (τ ) 2 + q 1 4E 0 + 4AE 0 I(q) (∂ τ θ q )(∂ τ θ -q ) dτ, (9)
where I(q) is the Fourier transform of the Coulomb potential, regularized 10 at very small q. ∂ τ θ q is in fact θq,j+1 -θq,j , and θq,j is the Fourier transform in space of θ i,j .

Rescaling the imaginary time and counting it in units of / √ 2E 0 J, and setting = 1 gives the effective action on a discrete space-time lattice

S ef f = N i=1 N ′ j=1 1 2 J 2E 0 θ i+1,j -θ i,j 2 + J 8E 0 q 1 1 + AI(q) N ′ j=1 θq,j+1 -θq,j θ-q,j+1 -θ-q,j , (10)
with N ′ = β √ 2E 0 J and θ the Fourier transform in space of θ i,j defined as

θq,j = 1 √ N N n=1 e 2πin θ n,j , (11)
for q multiples of 2π/N . Now, taking into account the periodic character of θ, we can take the effective action

S ef f = N i=1 N ′ j=1 J 2E 0 cos θ i+1,j -θ i,j + J 8E 0 q 1 1 + AI(q) N ′ j=1
θq,j+1 -θq,j θ-q,j+1 -θ-q,j .

Having derived an effective action, we first turn to the case where the periodic character of the phase can be neglected. The approximation is good for large J/E 0 . Cooper pairs feel very little Coulomb-blockade. The fugacity of phase slips goes essentially 3 as exp(-4 2J/E 0 ). Phase slips are neglected and a phase correlator derived. Owing to LRI, the system is no more a Luttinger liquid. Correlations ressemble some correlators appearing in one-dimensional quantum wires with low electron densities 9,10 , sometimes dubbed "Wigner crystals".

B. Spin waves

In this section, we allow only smooth variations of θ. It is then possible to expand the cosine term. This gives immediately a spin-wave Hamiltonian. We check here the known spin wave dispersion relations.

S ef f = ω q J 8E 0 q 2 + 1 1 + AI(q) ω 2 θ q,ω θ -q,-ω , (13)
with I(q) ≃ ln(q -1 ) for small q. Setting K = 1 2π 8E0 J , this results in

θ q,ω θ -q,-ω = πK q 2 + 1 1+AI(q) ω 2 , (14)
and thus

G(x, τ ) = θ(x, τ ) -θ(0, 0) 2 = 2πK 1 β ωn ∞ -∞ 1 -cos(ω n τ + qx) q 2 + ω 2 n 1+AI(q) dq √ 2π , (15)
where β = 1/(k B T ) and ω n are the bosonic Matsubara frequencies. For distances x and times τ not too long, so that the long range term does not influence the behavior, the correlations are those of a Luttinger liquid G(x, τ ) =

K ln √ x 2 +τ 2 α
, where α is a short range cutoff (corresponding to finite bandwidth), whereas for very large times 10 ,

G(x, τ ) = K √ 2A ln 3/2 x 2 + √ 2A τ √ ln τ 2 α , (16)
if τ is counted in units of / √ 2E 0 J. The correlation function of e iθ for spin-waves read g(τ ) ≡ e iθ(0,0) e -iθ(0,τ ) = exp -B ln(τ ) + 1 2 ln ln τ

3/2 , (17)
with B = K A 2 . We would like to calculate the full correlation function, not only in the spin wave approximation, but including phase-slips. This turns out to be more intricate, we will neglect all factors that vary slower than any power law, such as, for example, logarithmic corrections or terms as exp(-C 0 √ ln τ ) with C 0 a constant. Within this approximation, we have the spin-wave correlator

g(τ ) ≃ exp -B ln(τ ) 3/2 . ( 18
)
C. Phase slip Hamiltonian

Now, we take into account phase slips, that is vortices for the XY model in 1 space and 1 time dimension. The procedure to split the action into two parts, one due to spin-waves and the other to vortices is standard without LRI, however we recall it here with some details to point out the influence of LRI.

We use the model of Villain and transform the exponential of the cosine term in (12) as the sum over n integers but with a quadratic action 17 , involving θ and n. The space direction will be denoted as x and the imaginary time direction as y. The long range interaction occurs only in the x direction, for bonds directed in the y direction. Let us denote by

F (q) = 1 1 + AI(q) . ( 19
)
As a consequence, the Coulomb term, i.e. the second term on the r.h.s of Eq. ( 12), has to be modified from q F (q) N ′ j=1 θq,j+1 -θq,j θ-q,j+1 -θ-q,j to q F (q) j l θ l,j+1 -θ l,j -n y l,j e iql r θ r,j+1 -θ r,j -n y r,j e iqr , where n α r,j denotes the (2π times) integer attached to the bond which goes from site (r, j) to site (r, j+1), and pertaining to the y direction. The n y q,ky are simply the Fourier transforms of the n y r,m . The lattice spacing will be denoted a, and e x (e y ) will be the unit vector in the x (y) direction.

The part not involving the n k 's is of the type q F (q) ky θ q.ky θ -q,-ky so that finally the total action is of the form α=x,y k

F (q)δ α,y + δ α,x |iK α θ k -2πn α k | 2 + ǫ|θ k | 2 , (20)
where

K α = 2 sin(k.a α /2), (21)
and a x = ae x , a y = ae y where e x and e y are the unit lattice vectors. ǫ is an infinitesimally small positive quantity.

From now on, we will denote q = k x so that k has components (k x , k y ). The only difference with the normal XY model is that the coupling constant in the y direction has a dependence in k x for small k x (long wavelengths).

Setting

A x = J 2E0 and A y = A x F (q), S ef f = k α=x,y A α |K α θ k -2πn α k | 2 . ( 22
)
Diagonalizing this action yields the decomposition into a continuous field θk and a discrete field n k . The action must be of the form

S ef f = S θ + S V , (23)
S θ = k C(k) θ * k θk , (24)
S V = α,γ,k n k,α n * k,γ D α,γ (k), (25)
where α and γ take two values, x and y. The quantities D α,γ (k) and C(k) are functions of k. The linear transformation is given explicitly by

θk = θ k + 2iπ α A α K α n k,α ǫ + α A α K 2 α , (26)
C(k) = α=x,y F (k x )δ α,y + δ α,x , (27)
D α,γ (k) = 4π 2 A x A y A x K 2 x + A y K 2 y K α K γ (2δ α,γ -1), (28)
and δ α,γ = 1, if α = γ and zero otherwise.

Next, S V is rewritten in terms of vortices defined as

q k = i(K y n k,x -K x n k,y ), (29)
which yields the vortex Hamiltonian

H V = 4π 2 k =0 A x A y q k q -k A x K 2 x + A y K 2 y . (30)
The apparent divergence at k = 0 is in fact absent because of charge neutrality (q k=0 = i q i = 0). Going back to real space, denoting N tot = N N ′ , the total number of sites:

q ρ = 1 √ N tot k q k e ik.Rρ , (31)
enables to rewrite the vortex Hamiltonian as

H V = ρ,ρ ′ V ρ,ρ ′ q ρ q ρ ′ , (32)
with, taking into account values of A x , A y and F (q),

V ρ,ρ ′ = -4π 2 J 2E 0 1 N tot k 1 -cos k.(R ρ -R ρ ′ ) 1 [1 + AI(k x )]K 2 x + K 2 y . (33)
First, the self-interaction V ρ,ρ has to be withdrawn, but it is not so crucial because of charge neutrality ( ρ q ρ = 0). Second, the vortex Hamiltonian as it stands is incomplete. The vortex fugacity term, of the form ρ (ln y)q 2 ρ is missing, with ln y being the fugacity of a vortex. In fact, one should start with a generalized Villain model 18 , introduced by José, Kadanoff, Kirkpatrick and Nelson, hereafter denoted by JKKN. As noted in Ref. 18, such a term is naturally generated by real space (Migdal-Kadanoff) renormalization. It is therefore necessary to reinstall this term. Next, we look at the shape of the interaction between vortices.

D. Two regimes

Although it is known that, with LRI, the vortex-vortex coupling for large separations 15 (much larger than the lattice constant) goes only as in ln |R ρ -R ρ ′ |, we need to derive a quantitative expression. The behavior is dominated by the behavior of the integral for small |k|. For theses small wave-vectors, K x ≃ ak x and K y ≃ ak y , a being the lattice constant. We can take R ρ = 0 and R ρ ′ = ρ cos(ν)e x + ρ sin(ν)e y .

V ρ,ρ ′ ≃ -4π 2 I, I = J 2E 0 1 -cos k x ρ ′ cos(ν) + k y ρ ′ sin(ν) 1 + A I(k x ) k 2 x + k 2 y d 2 k. (34)
The cosine term will be close to 1 only for k x much lower than 1/[ρ ′ cos(ν)], and k y much lower than 1/[ρ ′ sin(ν)]. To extract the main behavior for large distances ρ ′ , we can replace it by the integral

I = 4 J 2E 0 km (ρ ′ cos(ν)) -1 km (ρ ′ sin(ν)) -1 1 1 + A ln(k -1 x ) k 2 x + k 2 y dk x dk y , (35)
where k m is a wave-vector much larger than 1 ρ ′ , but still sufficiently small that the expansion of K x and K y for small k holds. The integral reads

I = 4 J 2E 0 km 1/ ρ ′ cos(ν) 1 1 + A ln(k -1 x ) k x atan k m k x 1 + A ln(k -1 x ) -atan 1 ρ ′ sin(ν)k x 1 + A ln(k -1 x ) dk x . (36)
As ρ ′ gets large, for ν different from 0 and π/2, ρ ′ cos(ν)/ A ln(ρ ′ ) goes to infinity, so the first atan term in the bracket is close to π/2 whereas the second atan term remains much smaller than 1, provided tan(ν) is not too small. We can thus safely replace the term into brackets by π/2. The integral on k x is then readily performed and is in

A -1 1 + A ln(k -1
x ). The limit A = 0 is singular. It is better to discuss limiting cases. Case 1:

• ln ρ ′ cos(ν) ≪ A -1 ,
This corresponds to distances where the decrease of the correlations due to long range unscreened Coulomb interactions still has not taken place. We have then

V ρ,ρ ′ ∼ 2π J 2E 0 ln ρ ′ cos(ν) . ( 37
)
Case 2:

• ln ρ ′ cos(ν) ≫ A -1 ,
This corresponds to distances where the decrease of the correlations due to long range unscreened Coulomb interactions has taken place.

We have then

V ρ,ρ ′ ∼ 2π J 2E 0 ln ρ ′ cos(ν) /A. (38)
Thus, at very large distances, there is no longer sufficient interaction between vortices to hold them tied together, whatever the short range term is 15 . Mathematically the BKT transition is expected to disappear. However, we need to put this on a more quantitative basis. Furthermore, in real systems, the Coulomb interactions remain unscreened only up to a certain screening length ξ e , which can be large but is not infinite. Also, the magnitude of the LRI Coulomb term may be much smaller than the one of the local Coulomb term. In order to determine the phase diagram and the phase correlations, we need to implement a real space renormalization group scheme analog to the one developped in Refs. 2,18. Care has to be taken since LRI will eventually dominate under the RG group. The last term in Eq. ( 2) is more relevant than the second one. However, under the RG fow, it will take some time for the last term (LRI) to overcome the effect of the local one (i.e. the second term in Eq. ( 2)). The situation is reminiscent of LRI in one-dimensional quantum wires with very low electron densities; there are important differences however. There will be remnants of the BKT transition. To proceed, we perform a RG in real space closely following Kosterlitz'original paper 2 .

III. RENORMALIZATION GROUP IN REAL SPACE

A quick look at the situation is that, owing to the Coulomb potential, at very large distances R, the interaction between vortices goes like √ ln ρ (instead of ln ρ in the usual case) , and is too weak to overcome the entropy term for placing a free vortex, which goes as ln ρ. Nevertheless, the last terms in Eqs. (3,4) can be small compared to the second ones, even if they are less relevant in RG sense, their initial value can be much smaller. Thus, it makes sense to treat them as a small perturbation, at least initially. This leads us to investigate how the usual KT RG flow is modified.

A. Modified KT equations

Here, the partition function has a form which is slightly different from the usual KT flow. The gas is neutral, without electric field. n will be the number of "+" charges.

Z = n 1 n! 2 K 2n exp i,j J F(r i -r j )q i q j ln r i -r j τ , ( 39
)
where J is a coupling constant, J = 2π J 2E0 and F is a function that behaves as 1 for not too large r i -r j but decays as ln |r i -r j | -1/2 for very large |r i -r j |. A model F can be

F(r i -r j ) = 1 + A ln |r i -r j | -1 . ( 40
)
τ is the short time cutoff and K the fugacity.

In the case where the long range term is much smaller than the short range one (at the beginning of the renormalization procedure only), we can however follow the same steps as Kosterlitz. The details are shown in Appendix A. We express the RG equation, using the variables

X = J -2, (41)
Y = 4πKτ 2 2 , (42) z = ln τ. (43)
Now, with the long range interaction, this becomes

dX dz = -Y -B √ z, (44)
dY dz = -2XY, (45)
with

B = 4 3 √ A.
Note that B can be much smaller than one. Now, because B is non-zero, the line Y = 0 is no longer a line of (stable or unstable) fixed points.

B. Physical description of the flow

The flow is no longer autonomous, as soon as B = 0. We expect a deformation of the usual flow. In particular, the line Y = 0 which was a line of fixed points in regular KT flow is no longer a true line of fixed points.

If the LRI were absent, we could be either in the case where phase correlations are decaying exponentially (disordered phase) or in the "ordered" phase, where correlations decay as a power law. In the first case, adding the long range interactions will just make the phase correlation length smaller. The more interesting case is when we start from the ordered phase. We thus examine quantitatively how switching LRI drives us into the disordered phase and, in the next section, we shall examine what the phase correlations look like.

In the usual case, (without LRI), the flow is autonomous (i.e. does not depend explicitly on time) and trajectories are of the form X 2 -Y = X 2 0 -Y 0 ≡ C, where (X 0 , Y 0 ) is the starting point. Setting y ≡ √ Y , they are hyperbolae in the (X, y) plane. In the "ordered" phase, the flow on the invariant curves is: which integrates to

dY Y Y + (X 2 0 -Y 0 ) = -2 dz, (46)
2 √ C ln C Y (z) + 1 + C Y (z) -ln C Y 0 + 1 + C Y 0 = 2(z -z 0 ). ( 47
)
z 0 is the initial "time" which can be set to zero. Integration to obtain X is also straightforward. The important thing is that, as z becomes very large, Y , which is related to the fugacity, tends to zero exponentially, while X tends to a strictly positive value X f in , see Fig. 2.

In our case, (with LRI), we can see what happens by performing the following approximation:

-if B √ z is much smaller than Y , we neglect it, -when, it becomes of the order of Y , we do the opposite and neglect Y in the r. h. s. of the dX/dz expression, Eq. (44).

Thus, there are at least two stages. In the first stage, we have the usual KT flow whereas in the second stage the long range term comes into play. We shall see that there are three stages actually. We now describe the first stage. Starting from a definite value of Y , by the regular KT flow, Y will decrease. Suppose that Y has decreased by a substantial amount, such that Y will be much smaller than C. This enables to replace 1 + C Y (z) by C Y (z) . Inverting the relation yields

Y (z) = 4C D 2 exp - √ C(z -z 0 ) , (48)
with

D = C Y0 + 1 + C Y0 . ( C Y0
is not necessarily much larger than one). The physics is that C is proportional to the initial distance to the separatrix. It is proportional to (T c -T )/T c in the thermal X-Y model. Now, we have to stop when Y reaches the value B √ z -z 0 , because we are no longer allowed to neglect the term A √ z in Eq. ( 44). This happens when

4C D exp - √ C(z -z 0 ) = B √ z -z 0 . ( 49
)
To get the main feeling, we can show that the typical value of z for which this happens is z 1 , such that z 1 -z 0 ≃ 1 √ C . This corresponds to a typical imaginary time τ 1 , or corresponding lengthscale l 1 such that

τ 1 ≃ l 1 = 1 B 1/ √ C . ( 50
)
Neglecting ln(ln) contributions, Y a , the corresponding value of Y is

Y a = B 1 √ C ln( 1 B ). ( 51
)
Y a goes to zero as B goes to zero. Now, we will enter the second stage. We now describe this second stage, where the second term in Eq. ( 44) becomes larger than Y . Then, for z larger than z a , we integrate the RG equations, neglecting the Y term with respect to B √ z.

X = X a - 2B 3 z 3/2 , (52)
and, using the other flow equation to obtain Y ,

ln Y = ln Y a -2X a z + 8 15 B z 5/2 . ( 53
)
Then, X will decrease and Y will continue to decrease, till X passes the value 0. Y will be minimum for X = 0 and then will start to increase. In fact, it increases faster than exponentially, so it will quickly become larger than B √ z -z 0 , so that our approximation will fail. The value of z 2 for which this happens can be estimated, but anyway z 2 -z 0 ≪ z 1 -z 0 , corresponding to a length l 2 = e z2 ≪ l 1 .

Then, there is a third stage, where Y becomes larger than √ B(z -z 0 ) . Then, we can again neglect √ B(z -z 0 ) in front of Y . We are then in the third regime which has analogies with the KT flow but on the disordered part of the transition. This lasts till Y has reached an appreciable value such that one is no longer in the vicinity of the critical point. Let z 3 be the corresponding value of z, and l 3 = e z3-z2 . l 3 will be also of order

B -1 √ C . Since C is proportional to T c -T ,
where T c is the transition temperature without LRI, it is tempting to conclude that, since B goes as √ A, the correlation length with LRI ξ ′ behaves as

ξ ′ = D √ 2E 0 J exp ln A -1 /2 T c -T , (54)
with D an unimportant constant of order one. It has the same form as the KT correlation length in the disordered phase (without LRI) but ln(A -1 ) replaces the usual 2 constant b. Of course, ξ ′ tends to infinity as B tends to zero. The general shape of the "flow diagram" is best seen on Fig. 2. Renormalization has to be stopped when τ becomes equal to β . Equating ξ ′ to /(k B T ) yields a temperature T * given by

T * = (Dk B ) -1 2E o J exp - ln(A -1 ) 2 √ T c -T * , (55)
with T c the critical temperature in the absence of LRI. If experiments are performed at a temperature larger than T * , it will not be possible to see the influence of LRI.

Having the general shape of the "flow-diagram", we now turn to the determination of the phase correlation function.

IV. CORRELATION FUNCTIONS

This section is a little technical. Just as the Hamiltonian decouples into a vortex part, the phase correlation function also do. Estimation of the spin-wave part is straightforward, but the estimation of the vortex part is more intricate. In section A, we first determine the contribution of one vortex to the phase correlator. This involves a vortex influence function v(r). Comparison between the usual case and the case with LRI is emphasized.

Then, in section B, average must be carried out over the possible configurations of vortices, using the Boltzmann factor with the vortex Hamiltonian. A Coulomb gas formulation is necessary. Distinction must be done between smooth variations of the vortex influence function v(r), studied in section C, and abrupt variations (section D); there are thus two contributions. Summing both contributions gives the vortex part of the phase corrrelator. Adding the spin-wave part yields the total phase correlator in imaginary time. In section F, analytical continuation is performed to obtain it in real time. The phase correlator in real time is necessary for the calculation of the transport properties of JJ chains.

For the determination of the correlation function, we use the method of JKKN. In the vortex Hamitonian, the core energy has to added. Such a term naturally arises from their Migdal-Kadanoff renormalization group. Therefore, we must use the generalized Villain lodel. The corresponding vortex Hamiltonian is

H V = ρ,ρ ′ V ρ,ρ ′ q ρ q ρ ′ + (ln y) ρ q 2 ρ , (56)
with y the charge fugacity. The q ρ are integer charges living on the dual lattice of the original lattice. They represent vortices. The correlation function to be computed is e iθ(x,τ ) e -iθ(0,0) . Written in Coulomb gas language, this correlation function disentangles into a spin-wave part and a vortex part.

A. Determination of the vortex influence function

We closely follow the method 19,20 of Pelcovits. The first thing to do is to determine what is the influence of a vortex sitting at R on the vortex part of the correlation function g V (rr ′ ). The total correlation function is g tot (rr ′ ) = e iθ(r) e -iθ(r ′ ) .

(57)

Following JKKN, g tot (rr ′ ) is written as a product of nearest neighbor correlation functions, along a path P that connects r to r ′ on the real lattice. This is for the total correlation function (vortices plus spin-waves). After (Gaussian) integration of the spin-wave degrees of freedom,

g tot (r -r ′ ) = g SW (r -r ′ ) g V (r -r ′ ). ( 58
)
The vortex part, g V (rr ′ ) is given by

g V (r -r ′ ) = exp i R m(R) v(R) , (59)
with

v(R) = R ′ V (R -R ′ ) η l (R) -η r (R ′ ) , (60)
where R are all the positions on the dual lattice, the function η l (R) (η r (R)) is zero everywhere except on sites of the dual lattice that lie immediately on the left (right) of the path chosen joining r to r ′ . V (R -R ′ ) is simply V ρ,ρ ′ , see Eq. ( 34). The brackets mean thermodynamical average with respect to H V . This results directly from Gaussian integration, it can be understood as follows: the contribution of a vortex sitting at location R, and vorticity m(R) is exp i m(R) v(R) . The function v(R) depends on r and r ′ . It has also another physical interpretation. It is the phase generated at site R by a vortex-antivortex pair, where one vortex has m = 1 and is sitting on r and the antivortex (m = -1) is sitting on r ′ .

The determination of v(R) has been done in the literature without the long range Coulomb interaction (LRI). This function will be called u(R), and the potential energy between vortices will be denoted as V 0 ρ,ρ ′ . We explain how the LRI affect the results. Following JKKN, without LRI,

u(R) = R ′ G ′ (R -R ′ ) η l (R ′ -η r (R ′ ) , (61)
where For not too large distances, the interaction between vortices is still logarithmic and we just recall the known results for completeness. For simplicity, we take a path joining r and r ′ to be a straight line (see Fig. 3). Working in cartesian coordinates, and taking the origin at r, the x axis, along the line going from r to r ′ , the coordinates of R, r and r ′ will be (x, y), (0, 0) and (r ′ , 0) respectively. The current integration variable R ′ on the dual lattice immediately above the cut will be of coordinate (x ′ , y ′ ).

G ′ (R -R ′ ) = V 0 Rρ,R ′ ρ ′ -V 0 Rρ,Rρ ≃ V 0 ρ,ρ ′ -V 0 ρ,ρ . See Fig. 3.
u(R) = lim y ′ →0 r ′ 0 d dy ′ ln (x -x ′ ) 2 + (y -y ′ ) 2 dx ′ = atan(y/x) -atan y/(x -r ′ ) , (62)
which gives

u(R) = Φ(r -R) -Φ(r ′ -R), (63)
with Φ(Z) = atan(y z /x z ), if Z = (x z , y z ). Now, consider the LRI, at large distances

V ρ,ρ ′ = 2π J 2E 0 ln(ρ -ρ ′ )/A. ( 64
)
The function u(R) at large distances must be replaced now by v(R).

v(R) = lim y ′ →0 2π √ A J 2E 0 r ′ 0 d dy ′ ln (x -x ′ ) 2 + (y -y ′ ) 2 dx ′ . ( 65
)
Integrating by parts, to extract the main contribution at large distances, yields the result

v(R) = A -1/2 √ 2 J 2E 0 - 1 ln (r ′ -x) 2 + y 2 atan( r ′ -x y ) - 1 ln x 2 + y 2 atan( x y ) + A -1/2 √ 2 J 2E 0 F, (66)
with F a quantity which is smaller by a factor ln(x 2 +y 2 ) -1

than the main one. Thus, v(R) behaves 15 as ln(R)

-1 2
for large R.

B. High temperature expansion of the vortex part of the phase correlator

From Eq. ( 59), the average of exp i R m(R)v(R) has to be performed with the full vortex Hamiltonian (including core energy). The Boltzmann weight of a particular configuration m(R) is given by exp A m(R) with

A m(R) = R m 2 (R) ln y + R,R ′ m(R)m(R ′ ) V R-R ′ and thus g V (r -r ′ ) = {m(R)} exp A m(R) + i R m(R)v(R) {m(R)} exp A m(R) . ( 67
) A high temperature expansion of g V (r -r ′ ) is needed. Applying the Poisson formula f (m) = ∞ p=-∞ ∞ -∞ f (φ) e -2iπpφ dφ to the function f (m) = exp A m(R) + i R m(R)v(R)
and integrating over the continuous auxiliary field φ enables to write it as

g V (r -r ′ ) = Z ′ Z , (68)
with

Z ′ = {p(R)} exp - 1 2 R,R ′ v(R) + 2πp(R) v(R ′ ) + 2πp(R ′ ) m(0)m(R) cont , (69)
where the sum is over all configurations of p(R). m(0)m(R) cont is the correlation function calculated with weight A m(R) but as though the m were continous variables. Z is the same as Z ′ except that v(R) is absent. However, because a switch from the m(R) variables to the p(R) variables (through the Poisson formula) has been made, now a low temperature expansion is needed. A high T expansion in the m(R) variables transcribes into a low T expansion in the p(R) variables. In the absence of LRI, m(0)m(R) cont has been calculated in Refs. 21,22. Even with the LRI, the correlation m(0)m(R) cont remains short range and can be approximated by nearest-neighbor only. Substituting the exact m(0)m(R) cont by a nearest neighbor interaction results in

Z ′ = {p(R)} exp - 1 16π 2 2E 0 J R,R ′ v(R) + 2πp(R) -u(R ′ ) -2πp(R ′ ) 2 . ( 70
)
Following Refs. 19,20, a low temperature expansion is performed, keeping only the essential diagrams. To simplify further, only the non-overhanging self-avoiding-walks (SAW's) joining R to r ′ are kept which results in

g V (r -r ′ ) = n P n (r -r ′ )e -n/4 q 2E 0 J exp 1 8π 2 2E 0 J ′ |∇ R v(R)| 2 d 2 R , (71)
where ′ means integration over all the plane except on the SAW. P n (rr ′ ) is the number of non-overhanging SAW's of n steps joining r and r ′ .

C. Smooth variations of v(R)

The calculation of the correlation function splits into the integration on the smooth variations of the field v(R) and on the abrupt discontinuities caused by the different diagrams occuring in the high T expansion of Z ′ . In our case, the situation looks a little more complicated because v(R) does not satisfy the Laplace equation. However, this causes no trouble because the volume contributions coming from ∆v, where ∆ is the Laplacian, will turn out to be very small. The main contributions come from the regions near the path P going from r to r ′ .

More precisely, the weight

W ≡ exp - 1 8π 2 2E 0 J ′ ∇ R v(R) 2 d 2 R , (72)
is to be evaluated, where the prime means integrations everywhere in the plane except on the diagrams. In the usual case, integration is made by parts

′ ∇ R u(R) 2 d 2 R = T 1 + T 2 + T 3 , (73)
T 1 ≡ - ′ S u ∆u d 2 R, (74)
T 2 ≡ Γout u ∇u .dσ, (75)
T 3 ≡ Γin u ∇u .dσ. (76)
S denotes the shaded aera in Fig. 4, Γ in is a contour which encloses the cut (from r to r ′ ) at a small distance b 0 and Γ out is a circle of radius R where R will eventually be allowed to go to ∞. dσ is the unit normal to the graph. On Fig. 4, it is simply the normal to the line joining r to r ′ . The first term T 1 vanishes since ∆u = 0. The second term T 2 will tend to zero as R goes to infinity. The last term T 3 has to be calculated using Eq. ( 61). Care has to be taken that G ′ (R -R ′ ) behaves only as ln(R -R ′ ) for sufficiently large distances and this form of G ′ should not be used for small distances. With this in mind, taking r = 0 and thus r ′ = |r ′ -r|,

T 3 = 2 r ′ 0 atan( y ′ x ′ ) -atan( y ′ r ′ -x ′ ) d dy ′ atan( y ′ x ′ ) -atan( y ′ r ′ -x ′ ) dx ′ . ( 77
)
The factor 2 comes from the fact that integration above and below the cut r, r ′ has to be performed (see Fig. 4). We take the integration path a finite amount b 0 above (or below) the cut to avoid spurious divergencies (see above). Then b 0 will be unimportant. The angle difference in the difference of the two atan functions Eq. (77) will correspond to the change in angle, taken to be π. Dropping the spurious ln(b 2 0 ) terms yields [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF] exactly T 3 = 4π ln |r ′ -r| . Now, we turn to the case of LRI. The analog of the terms T i , i = 1 to 3 will be denoted as T ′ i and u has to be replaced by v. Calculations are tedious and relegated to the appendices. It turns out that T ′ 1 and T ′ 2 are completely negligible whereas T ′ 3 = B 1 ln ln |rr ′ | , with B 1 a constant. Even after taking the exponential of this will only lead to logarithmic factors for the correlation functions and we thus omit this factor.

D. Abrupt variations of v(R), graph contributions

Here, the part due to abrupt variations of the field v(R), is the same as without LRI sum on SAW joining r and r ′ and the result is analogous to the usual case, namely, the contribution, called g V SAW reads:

g V SAW ≃ exp - |r-r ′ |ωp ξ |r -r ′ |ω p , (78)
ω p is the plasma frequency √ 8E 0 J/ . Being the result of a high-T expansion, this expression corresponds to short correlation length ξ. It is probably not valid for very large ξ. In the limit A → 0, the usual divergence close to the BKT transition must be retrieved.

E. Final form

The spin-wave part, the vortex part coming from smooth variations of v(R) and the vortex part coming from abrupt variations of v(R) just need to be multiplied. Thus we have, with r -

r ′ = |r -r ′ |, g(r -r ′ ) = exp -B ln |(r -r ′ )ω p | 3/2 exp(B 1 ln ln |(r -r ′ )ω p |) exp(- |(r-r ′ )ωp ξ ) (r -r ′ )ω p , (79)
with B 1 a constant. Contrary to the usual case, the part of the vortex contribution corresponding to smooth variations of v(R) no longer cancels the spin-wave contribution. Essentially

g(r -r ′ ) = exp - |r -r ′ |ω p ξ exp -B ln |(r -r ′ )ω p | 3/2 |(r -r ′ )ω p | -1/2 , (80)
where we systematically dropped all contributions going slower than any power law. We see that, contrary to what happens in the usual case, the total correlator is not of the Ornstein-Zernike type. However, in the limit B → 0, the Ornstein-Zernike form is retrieved.

F. Analytic continuation

It is now necessary to obtain the phase correlator in real time, and then in real frequency. We need to perform an analytical continuation τ -→ i t. We do not attempt to get the short time behavior. This would correspond to high voltages. At long times,

p(t) ∝ exp(-|ω p t|ξ -1 ) |ω p t| -1/2 exp -B ln |ω p t| 3/2 , (81)
where ∝ means "proportional to". Normalization of p would require more information on the shape of p(t) for shorter times. This yields in real frequencies Ω, for small Ω's (with respect to ω p )

P (Ω) ∝ Θ(Ω -ω p ξ -1 ) exp -B ln |Ω/ω p -ξ -1 | 3/2 (Ω/ω p -ξ -1 ) -1/2 , ( 82
)
where Θ is the Heaviside function. For larger Ω, this form is not valid. One would need the short time behavior of p(t), which we do not have explicitly. However, on those short timescales, the long range part of the interaction will play a minor role. Therefore, in order to perform calculations, a reasonable assumption is that, as without LRI, P (Ω) is negligible 14 for Ω larger than κω p , with κ = 2. Enforcing the normalization ∞ 0 P (Ω) dΩ = 1, we thus have the approximate form

P (Ω) = p 0 ω -1 p Θ(Ω/ω p -ξ -1 ) exp -B ln |Ω/ω p -ξ -1 | 3/2 (Ω/ω p -ξ -1 ) -1/2 , (83)
with dx.

p 0 = 2L(2 3/2 B, κ -ξ -1 ) -1 ,

V. PHYSICAL APPLICATIONS

Having obtained the phase correlator in real time, we then use it for the calculation of the Andreev current. A normal metal is placed on the left of the Josephson junction chain and biased by a voltage V , below the gap |∆|, see Fig. 1.

A. Andreev average current

The average Andreev current flowing from the metal into the superconductor can be readily expressed as a function of the distribution of phase modes P (Ω). The case P (Ω) = δ(Ω) corresponds to a usual BCS superconductor, and BTK results [START_REF] Blonder | [END_REF] are retrieved. The case where P (Ω) is given by the Luttinger liquid theory or usual 2d XY model has been studied in Refs. 14,26. The contact between the metal and the JJ chain can range from tunnel to perfect contact. Following the Keldysh formalism, and performing a nonperturbative expansion, to all orders in the tunneling amplitude, one can describe all situations 25 . Without phase fluctuations, the procedure has been shown to be exactly equivalent to solving the Bogolubov-de Gennes (BdG) equations. With phase fluctuations, it is only an approximation. The average current 26 reads:

I = 16e|Γ| 4 h W 4 |∆| 2 2 eV 0 P (Ω) dΩ eV -Ω 2 0 D -1 (ω) dω, (84)
with

D(ω) = |∆| 2 2 -ω 2 1 + |Γ| 4 W 4 2 + 4 |Γ| 4 W 4 ω 2 , (85)
W being the bandwidth and Γ a hopping parameter. For tunnel contacts, |Γ| ≪ W . A perfect contact corresponds to Γ = W exactly. See Ref. 25 for details. We would like to compare two cases, with and without LRI. For JJ chains, LRI are usually not taken into account. One exception are the models of JJ chain considered in Ref. 27, where long range interactions in the capacitance matrix were examined. These models yield physical results which are esentially the same as the usual XY model, except in the case where capacitance of the islands to the ground, C -1 g is strictly zero. Only in this case does the BKT transition disappear and the chain is always in the disordered phase, because of the proliferation of unbound vortices.

In most experimental situations, it does not seem to be realistic to assume that C -1 g is really completely negligible with respect to nearest-neighbor islands capacitive coupling (or even with respect to longer range inter-island capacitive couplings). One exception may be the case of the experiments of Ref. 6. By contrast, the model studied in this paper always keeps a non-zero capacitance to the ground and still always gives unbound vortices. Now, we would like to compare the results of our calculation with a model commonly used in the literature. In this latter model, a nearest-neighbor interaction bewteen island charges is included, and there is no LRI. We discuss qualitatively this case and compare the predictions from such a model with our case. A term of the form i q i q i+1 /(2C 1 ) must be added to the Hamiltonian, see Eq.( 3). Qualitatively, this will add some charge fluctuations and should decrease the phase correlations. In a more quantitative way, this will change for example the correlator of θ in Fourier space Eq. ( 14) to

θ q,ω θ -q,-ω = πK q 2 1 + C C1 cos(qa) + ω 2 , ( 86
)
where a is the spacing between Josephson junctions. Calculations are straightforward and the main changes with respect to the case with just a local Coulomb term are the following. In the quasi long range phase, for frequencies much smaller than ω p , the plasma frequency is now ω ′ p = ω p 1 + C C1 , and thus the constant K has to be changed accordingly to

K ′ = K 1 + C
C1 . This no longer holds for frequencies of order ω p . The critical point will be shifted. Defining t ≡ ω p /J, the usual BKT transition which occurred at t = t c = π/2 will now take place at

t ′ c = π 2 q 1+ C C 1 .
In the disordered phase, but still no too far from criticality, the correlation length will be depressed. However, the phase correlator will retain approximately its Ornstein-Zernike form, Eq. ( 78), with different values of ω p and ξ. The important thing is that the pre-exponential factor is still in 1/ |rr ′ |.

As an example, we take a perfectly transmitting interface and a situation where the Josephson chain would be insulating (disordered phase) even in absence of LRI. The presence of LRI will modify the correlation length ξ but also alter the form of the correlator P (Ω). Then, D(ω) is a constant and the current assumes the form

I = p 0 eω p π X 0 exp -B| ln x| 3/2 √ x (X -x) dx, (87)
with

X = 2eV ω p -ξ -1 . (88)
I remains zero for voltages smaller than V 0 = ωp 2e ξ -1 , just as without LRI. The new thing is that, for voltages slightly larger than V 0 , I no longer behaves as (V -V 0 ) 3/2 . We introduce the dimensionless quantity l -1 = exp(-B -2/3 ) and the corresponding voltage V 1 = ωp 2e (l -1 + ξ -1 ). V 1 is a crossover value of the voltage, not a sharp threshold as V 0 . For voltages V , such that V 0 < V ≪ V 1 , then, X ≪ l -1 and I behaves as

I = C 2 V -V 0 V 0 exp -B ln V -V 0 V 0 ξ -1 3/2 , ( 89
)
with C 2 a constant. Thus, the current starts to increase faster than any power law of (V -V 0 ). If, on the other hand, V is larger than V 1 , then, I increases again as (V -V 0 ) 3/2 . For even larger voltages, properties of the phase correlator at short time scales are needed. We expect that, the LRI will not be relevant on these scales and that the linear behavior of the I -V curves in the usual case will be retrieved. This is illustrated in Fig. 5, where we plot I in units of eω p /π versus eV /(2 ω p ), for t = 2π. Since our calculations describe only moderate LRI, for the LRI parameter A, we take A = 0.2 which results in B = 0.316. For the inverse correlation "length" ξ -1 , we take 2 it to be exp -1.5/ √ t -t c . This serves just as a rough estimate since the exponential behavior of the correlation length applies only (above and) close to the usual BKT transition point. It leads to ξ -1 /4 ≃ 0.125, to be compared to l -1 /4 ≃ 0.029. We see that the Andreev current is strongly depressed by LRI. We also compare with the model without LRI but with nearest neighbor Coulomb interactions, mentioned above. We have taken C C1 = 0.25. This causes the voltage threshold V 0 to be increased to a slightly higher value V ′ 0 , but the current still starts as (V -V ′ 0 ) 3/2 ; see Fig. (5).

To be more quantitative, we make a logarithmic plot, on Fig. 6, of I versus (V -V 0 /V 0 ) ξ -1 . Whereas the curve without LRI, for A ≃ 0 has a slope 3/2, this is not the case with LRI. The solid curve catches only with the curve without LRI, for V -V 0 of the order of 0.1 to 1, which is slightly larger than of order l -1 /4. For the model with nearest-neighbor interaction, we plot instead I versus (V -V ′ 0 )/V ′ 0 ξ ′ -1 . The curve (dash-dotted line) is still very close to a straight line of slope 3/2.

In the case where the Josephson chain would be in the (quasi) ordered phase, t < t c , without LRI, LRI impose a correlation length ξ 1 of order

ξ 1 = A - 1 2 √ tc -t , (90)
and thus a threshold voltage ( ω p /2e) ξ -1 1 .

eV/2 h ---ω p 〈 I 〉 0.0 0.2 0.4 0.6 1.

0.

FIG. 5: Andreev current I , [ in units of eωp/π ], as a function of applied voltage V , [ in units of 2 ωp/e ], for a perfectly transmitting interface and ωp/J = 6.28, with LRI and A = 0.2 (solid line), without LRI (A = 0) but with a nearest-neighbor interaction (dash-dotted line), with only on-site interaction (dashed line).

[

(V-V 0 ) /V 0 ) ] ξ -1 I / / / / 10 - 4 10 -3 10 -1 1. 10 -8 1. 10 -2
FIG. 6: log -log plot of Andreev current I , as a function of h (V -V0)/V0 i ξ -1 , (see text), same legend and values as Fig. 5, except for the case with a nearest-neigbor interaction (dash-dotted line). In this case the threshold value V ′ 0 is different from the threshold V0 for the two other cases and we plot h

(V -V ′ 0 )/V ′ 0 i ξ ′ -1 instead on abscissa.

VI. CONCLUSION

A study of the influence of unscreened long range Coulomb interaction was carried out for one-dimensional array of Josephson junctions. This can also have implications for one-dimensional thin superconducting wires. We focussed on the case where the initial term describing LRI was much smaller than Josephson coupling and local Coulomb blockade term. Deriving a Hamiltonian describing the interaction of phase slips, we examined how the KT flow is modified by LRI.

Implications on Andreev reflexion between a normal metal and a superconductor was examined and LRI introduce a new voltage scale. The analytical form of the I -V characteristics was derived for any value of the transparency of the normal-superconductor interface, (from tunnel contact to perfect interface). If we start from a situation where the charging energy dominates the Josephson coupling, corresponding to the JJ chain being insulating, the introduction of LRI changes the behaviour above the threshold voltage V c from a (V -V c ) 3/2 behavior to a sharper increase, faster than any power-law of (V -V c ). This lasts up to a crossover voltage V 1 , whose value depends on the strength of the long range Coulomb interactions. If, on the contrary, we start from a situation where the Josephson coupling dominates without LRI (quasi ordered phase), then LRI will always generate a threshold voltage V 0 below which no Andreev current can flow. This threshold V 0 goes as A 1 √ tc -t where A is proportional to the strength of the LRI and t -t c is the distance to the critical point in the absence of LRI, with t = ω p /J and t c = π/2. With LRI, above V 0 , the current starts to increase faster than any power law. The form of the phase correlator obtained in this paper could be used for calculation of other physical quantities relevant to JJ chains or to one-dimensional superconducting thin quantum wires.

Such calculations are not directly relevant for high-T c superconductors 13 , because of their layered structure and their strong anisotropy in the direction transverse to the planes, whereas the former calculation remained striclty one-dimensional.

For designed JJ chains, the inverse capacitance is generally cut exponentially with a screening length of the order of λ s ∼ a C/C g , with usually λ s ≃ 30a, (a is the spacing between neighboring junctions), [START_REF]Single Charge Tunneling[END_REF][START_REF] Odintsov | [END_REF] . The number N of junctions 34 can be up to 255, so that λ s ≪ L ≡ N a. Thus, the model considered here should not be construed as a model for fabricated JJ chains but rather an attempt to include LRI, which do exist in isolated thin wires. However, we can use the above calculations, replacing the badly screened Coulomb potential by the exponentially screened potential U ij , see Eq. ( 2). This leads qualitatively to an interaction between vortices V ′ ρ,ρ ′ which has the following properties. For distances ρ much smaller than λa, V ρ,ρ ′ = D J 2E0 λ -1/2 (a/ρ) ln ρ, with D a constant of the order of ten. The strength of the interaction between vortices decays very rapidly, (much faster than for the badly screened Coulomb interaction). This lasts till ρ becomes of the order of λa, where V ′ ρ,ρ ′ stops to decrease and behaves as V ρ,ρ ′ = D J 2E0 λ -3/2 Lnρ, for large distances, i.e. ρ ≫ aλ, Theoretically, the RG group is started at coarse graining scales smaller than λa and as the RG flow proceeds, the interaction between vortices gets smaller. This may tend to drive towards the disordered region. This is not unexpected, given the large amount of Coulomb repulsion provided by the charge coupling of the islands which are less than λa apart. When RG proceeds further and the coarse-graining length becomes larger than λa, then, the interaction V ρ,ρ ′ bewteen vortices ceases to decay and becomes logarithmic with distance. According to the value of the coupling reached at that moment, we can flow either to the disorded region or to a quasi-ordered fixed point. The screening of the Coulomb interaction renders quasi long range-ordering more favorable. In particular, contrary to the badly screened interaction case, phase correlations do not always decay exponentially at the end.

Let us suppose first that J is sufficiently large so that we end in the (quasi)-ordered phase, P (Ω) will decay as a power law for small Ω, with some exponent α, for Ω smaller than ω p /λ. For Ω larger than ω p /λ, the decay of P (Ω) corresponds to a larger effective vortex interaction, and phase correlations are enhanced, leading to a larger P (Ω). As we go to large distances (or longer times or smaller voltages), vortex interactions get reduced until the distance λa is reached. This is due large Coulomb repulsion at shorter lengthscales. This has experimental consequences on the shape of I vs. V curves. Frequencies smaller than ω p /λ correspond to voltages smaller than ω p /e. For these voltages, I -V curves will start as power law with some exponent, normally as α with α = 1 + 1 π E ′ 0 J ′ ), where E ′ 0 and J ′ are renormalized values of E 0 and J. When V becomes larger than ω p /e, the current is enhanced and increases faster than V α . This remains true for voltages still smaller than ω p /e. Now, let us suppose that J is not so large and one ends in the disordered phase. There is a threshold voltage but beyond V c , the current starts to increase as (V -V c ) 3/2 in the usual theory. Inclusion of screened interaction leads to an increase which is faster than (V -V c ) 3/2 , which could also be tested experimentally.

The model with badly screened interactions may also have applications for very long (but not unrealistically long) 1d superconductor thin wires. More precisely in thin wires (typically 3 to 5 nm diameter), the interaction between phase slips has the same structure as the Hamiltonian (1), but the physical origin is somewhat different 28 . First, in concrete situations, it is not the two-dimensional BTK transition which is seen but rather a one-dimensional Schmid transition 31 . We briefly recall the reasons below.

In experiments with very homogenous thin wires in the difusive limit, of diameter 3 to 10 nm, such as the ones used in Refs. 29,30 the maximum length L is about 1µm. Operating temperatures are between 100mK and a few Kelvin. If the length were infinite, renormalization had to be stopped when coarse graining length τ on the imaginary time axis reaches β . However, in practise, two-dimensional renormalization has to be stopped when it reaches L/v MS where v MS is the Mooj-Schön velocity 7 . v MS ∼ ω p d/C L is equivalent of the spin-wave velocity, where C L is the capacity per unit length and d the diameter of the wire. v MS is of the order 10 5 m.s -1 , so that for a 1µm long wire, L/v s is of the order 10 -11 s and for temperatures lower than 0.75K, renormalization has to be stopped first on the x axis, when the coarse graining length reaches the size of the wire. We end up with charges in 1 dimension (as for a single Josephson junction) which interact logarithmically as I ren Ln|r i -r j |, with I ren a renormalized coupling. Then, renormalization proceeds solely on the imaginary time direction. Depending on the value of I ren , fugacity of vortices renormalizes to zero or not, resulting in the usual Schmid transition. The spatial structure on the x axis has been completely lost and the long-range interaction do not play any role, except that they changed slightly the values of I ren . Thus, in this case, LRI do not seem to play an important role. Thus, in order to see the influence of LRI, longer wires are needed.

For longer wires (but not totally unrealistic, a few microns long), we can expect temperatures larger than v s /(Lk B ), so that renormalization has to be cut on the imaginary time axis first. Then, LRI will drive the resistance R to a finite value. If R versus T curves are plotted, no curve should plunge towards zero resistance, unless T becomes of the order of v s /(Lk B ), which is about 0.5 K for L = 10µm.

There are difficulties to see LRI. One could be that stray capacities (unavoidable in experiments) can also cut the long-range term in the capacitance for small q (in real space) so that C -1 (k) is no longer of the form ln k for small k. This depends on the charateristics of the wires. Another problem is that the action for a single phase slip, which controls its bare fugacity, is strongly dependent on temperature and on the diameter of the wire. This dependence seems to be more important than the interaction bewteen phase-slips in determining the shape of the curves showing the resistance R versus T , for different wire thicknesses 28 .

Long range interactions are originally not screened in isolated one-dimensional wires. However, the presence of an environment could alter the screening. Situations where a plate or a wire is placed in the vicinity of a superconducting thin wire are now experimentally 35 and theoretically studied 8,36 . This aims at providing mechanisms of dissipation so as to damp the phase fluctuations. Though we did not consider in this paper coupling to an external environment, our calculation may be viewed as a first step towards understanding the role of LRI, including, both spin waves and phase slips (solitons).

FIG. 2 :

 2 FIG.2: Modified KT flow in the presence of long-range interactions. Only a point which would correspond to the "ordered" phase of the usual KT diagram is shown. Here, because of long-range interactions, we always end in a completely disordered phase.

FIG. 3 :

 3 FIG.3: Determination of influence of vortex sitting at R on the correlation function e iθ(r) e -iθ(r ′ ) .

FIG. 4 :

 4 FIG.4: Determination of influence of the smooth variations of the field u(R) on the vortex part of the correlation function e iθ(r) e -iθ(r ′ ) .

  and the function L is defined by L(b, y) = y 0 exp -b ln |x| 3/2

Acknowledgments

We thank A. Crépieux and J. Rech for reading of the manuscript.

APPENDIX A

In this appendix, we derive the RG flow in the presence of long-range interactions. We define g as g(R i,j ) ≡ J F(r ir j ), (A1) with

Integration has to be done over r j in the whole plane except the two circles of radii τ around r k and r l . Setting r = (r l + r k )/2 and denoting by K the vortex fugacity, the partition function Z reads

where D is the whole plane except aera inside the circle of radii τ around the r i , i = 1, n. Let us denote by

where D is the whole plane except the two circles of radii τ around r k and r l . Given the value of F(r), see Eq. (40), we have

for A ≪ 1. Denote by

where S is the aera. Collecting all the terms,

Using e x = 1 + x, for small x,

When changing τ into τ + dτ , this leads only to the renormalization of the fugacity K into K 1 -J dτ τ replacing F(r i -r j ) by its initial value, which is one. After changing (k, l) → (i, j) in the second summation, the part in the exponential reads i,j

If F(R i,j ) were a constant, we would retrieve the usual KT flow. However, since we are looking at larger and larger lengthscales, the effective coupling will not only decrease because of RG but also because of its √ ln R dependence. The smallest R contributing to Z is τ . Thus, adding a term proportional to - 4 3 A ln(τ ) will take this into account. Now, the partition function has to be rearranged so as to look as its original form, with renormalized parameters and an overall multiplicative factor, setting z = ln τ , we obtain

Switching to the usual variables

and linearizing around X = 0, the flow reads

Here, we give details on the contributions to the vortex part of the correlation function g V (rr ′ ) coming from smooth variations of v(r). The first term to evaluate is

where the double integral extends over the whole shaded surface in Fig. 4. ∆ is the Laplacian operator (with respect to variable R). In this calculation, r ′ and r are fixed, and, without restriction, r is taken to be zero for simplicity. If R is the modulus of R, since v decreases as ln(R) -1/2 for large R, it turns out that ∆v behaves, as R -2 (ln R) -3/2 . Thus, v ∆v behaves as R -2 (ln R) -2 . When surface integration is performed, the integral tends to a constant, for r ′ → ∞ because of the (ln R) -2 factor. When exponentiated this will give a constant. Since, we shall neglect all terms that depend on r ′ , weaker than any power law, we shall neglect this factor.

The second term is a line integral over the outer contour Γ out on Fig. 4,

where dσ is the normal to the contour Γ out which is here chosen as a circle of radius R. For all R's belonging to Γ out , R are equal to R and v are of order (ln R) -1/2 , the leading terms in (∇v) are of order R -1 (ln R) -1/2 . Thus T 2 is at most of order (ln R) -1 and vanishes as R goes to infinity. The third term to evaluate is a line integral over the inner contour Γ in . on Fig. 4.

∇v).dσ. (B3)

A calculation completely analogous to the case without LRI shows that the leading terms are of order

x ′ 2 +y ′ 2 ln(x ′ 2 + y ′ 2 ) -1

dx ′ + x ′ ↔ (r ′ -x ′ ) , with r ′ the modulus of rr ′ . This behaves, for large r ′ as ln(ln r ′ ).