Influence of vortices and phase fluctuations on thermoelectric transport properties of superconductors in a magnetic field

Influence of vortices and phase fluctuations on thermoelectric transport properties of superconductors in a magnetic field

Andreas Andersson    Jack Lidmar Theoretical Physics, Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden
July 12, 2019

We study heat transport and thermoelectric effects in two-dimensional superconductors in a magnetic field. These are modeled as granular Josephson-junction arrays, forming either regular or random lattices. We employ two different models for the dynamics, relaxational model-A dynamics or resistively and capacitively shunted Josephson junction (RCSJ) dynamics. We derive expressions for the heat current in these models, which are then used in numerical simulations to calculate the heat conductivity and the Nernst coefficient for different temperatures and magnetic fields. At low temperatures and zero magnetic field the heat conductivity in the RCSJ model is calculated analytically from a spin wave approximation, and is seen to have an anomalous logarithmic dependence on the system size, and also to diverge in the completely overdamped limit . From our simulations we find at low magnetic fields that the Nernst signal displays a characteristic “tilted hill” profile similar to experiments and a nonmonotonic temperature dependence of the heat conductivity. We also investigate the effects of granularity and randomness, which become important for higher magnetic fields. In this regime geometric frustration strongly influences the results in both regular and random systems and leads to highly nontrivial magnetic field dependencies of the studied transport coefficients.

74.81.-g, 74.25.F-, 74.25.fg, 74.25.Uv, 74.78.-w

I Introduction

Thermoelectric effects in superconductors are of considerable interest, since they provide an important probe of fluctuations and correlations in these materials. Such effects have gained an increasing amount of attention since the recent measurements of the Nernst effect in the pseudogap phase of underdoped high- cuprates, where a remarkably large effect was observed Xu2000 (); *WangOng2006. The Nernst effect has since then been measured in many other materials, e.g., in superconducting thin films Pourret2006 (); *PourretBehnia_NJP2009 and in the iron-pnictides Zhu2008 (). The Nernst effect is usually very small for conventional metals, making it a particularly sensitive probe of superconducting fluctuations. Theoretical and numerical studies have described the large Nernst effect either in terms of superconducting fluctuations of Gaussian nature Ussishkin (); MukerjeeHuse (); Serbyn2009 (); Michaeli2009 (); *Michaeli2009a, or as fluctuations of the phase of the order parameter, i.e., as a result of vortices WangOng2006 (); Podolsky (); Raghu2008 (). Other explanations of nonsuperconducting origin have also been put forward, e.g., proximity to a quantum critical point Bhaseen2007 (); *Bhaseen2009; Sachdev (), quasiparticles Oganesyan2004 (); Zhang2010 (), stripe order Hackl2010 (), etc. Here we will focus exclusively on the effects of phase fluctuations and vortices on heat and charge transport. Phase fluctuations were early on proposed to play a key role in the pseudogap phase of underdoped cuprates Emery1995 (). In quasi-two-dimensional superconducting films and Josephson-junction arrays they are known to be the dominant fluctuations Berezinskii1 (); *Berezinskii2; *KT.

In this paper we study the heat transport and the thermoelectric response in two-dimensional granular superconductors or Josephson- junction arrays, using two widely employed models for the dynamics, (i) relaxational model-A dynamics, described by a Langevin equation, or (ii) resistively and capacitively shunted Josephson junction (RCSJ) dynamics. These are well suited for numerical simulation studies, and have been used extensively to calculate electric and magnetic static and dynamic properties, and to study vortex dynamics under influence of electric currents MonTeitel (); ChungLeeStroud (); Kim (); Marconi (). The calculation of thermoelectric properties is, however, less straightforward. Previous simulation studies have used time-dependent Ginzburg-Landau theory MukerjeeHuse (), phase-only XY models Podolsky () with Langevin dynamics, or vortex dynamics Raghu2008 (). They have been limited to rather narrow parameter regimes with a focus on explaining the large Nernst effect seen in the pseudogap phase of underdoped cuprates. Here we present a comprehensive study of heat conductivity, thermoelectric effects, and resistivity for a broad range of parameters. We also investigate the effects of a granular structure. The models are defined on a discrete lattice and can be experimentally realized in granular superconductors and artificially fabricated Josephson-junction arrays. At low magnetic fields discreteness effects become less important so that our results in this regime are relevant also for homogenous bulk superconductors. On the other hand, in granular superconductors transport properties are strongly affected by discreteness and geometric frustration effects Beloborodov2007 (). This is particularly true at high magnetic fields when the inter-vortex distance becomes comparable to the granularity. This leads to a rich structure in, e.g., the Nernst signal as the magnetic field is varied, with anomalous sign changes occurring in the vicinity of special commensurate fillings AnderssonLidmar ().

To start with, let us first recall that the heat current density and electric current density are related to the thermodynamic forces, the electric field and the temperature gradient , by the standard phenomenological linear relations


where the thermoelectric and electrothermal tensors obey the Onsager relation . The Nernst coefficient is defined as the off-diagonal response of the electric field to an applied temperature gradient in a transverse magnetic field ,


where is the so called Nernst signal. Both the Nernst effect and the heat conductivity are measured under the condition of zero electric current, such that


In a system with no Hall effect (), which is the case for the models studied below, Eq. (3) reduces to .

In a phase-fluctuating superconductor, mobile vortices, either thermally excited or induced by an applied magnetic field, may significantly affect transport properties. An applied electric current will exert a perpendicular force on the vortices and their motion will generate a transverse electric field parallel to , thus causing a finite flux-flow resistivity. Vortex motion can also be caused by a temperature gradient, thereby generating an electric field perpendicular to both the magnetic field and the temperature gradient. The Nernst coefficient defined in Eq. (2) can be seen simply as the diagonal response of the vortex velocity to the temperature gradient. For this reason it is plausible that a large Nernst signal is expected in the vortex liquid phase.

Via an Onsager relation a heat current can then be generated from an applied electric current. The vortices therefore also contribute to the heat conductivity, although other heat carriers, e.g., phonons or quasiparticles, probably dominate. From the vortex point of view it is natural to consider the applied current as the driving force and the electric field, which is proportional (but perpendicular) to the vortex current, as the response. This yields an alternative, but equivalent formulation of the linear relations in Eq. (1)


where . This is the approach employed in our simulations. Instead of calculating through other transport coefficients, i.e., using Eq. (3), we obtain the Nernst signal directly as for . The longitudinal heat conductivity is just the diagonal components of the tensor in Eq. (5), and is the resistivity.

The picture described above applies when the vortices are free to move in response to the driving forces. Pinning of vortices to material defects, grain boundaries, etc., can lead to dramatic changes of the transport properties.

The paper is organized as follows. In Sec. II we introduce the models and their dynamics, Langevin or RCSJ, on general two-dimensional (2D) lattices. In Sec. III we derive an expression for the heat current for the models studied. This has over the years proven to be a subtle task, especially in the presence of a magnetic field. We present two separate routes to finding the explicit form of the heat current for Langevin and RCSJ dynamics. In addition we show, using a functional integral approach, that the Nernst signal indeed can be calculated from a Kubo formula involving the heat current. Section IV discusses the thermal conductivity at zero magnetic field in the low temperature limit, where a spin wave approximation is applicable. Our analytic calculations reveal a logarithmic system size dependence of in this regime for the RCSJ case. In Sec. V we give some technical details of our numerical simulations, and in the last part, Sec. VI, the results of our numerical simulations on square and irregular lattices are presented. We consider the case of zero and weak magnetic fields as well as the high field limit, where the transport coefficients are severely affected by geometric frustration. In the weak magnetic field limit the results are discussed in relation to previous theoretical works and experiments.

Ii Models

We model a 2D granular superconductor (of size ) as an array of superconducting grains connected by Josephson junctions. These grains may or may not be ordered and connected in a regular fashion. The supercurrent flowing between two superconducting grains and is given by the Josephson current-phase relation


where is the critical current of the junction, is the superconducting flux quantum, and is the superconducting phase of grain . The grains are assumed to be small enough ( the coherence length) that the phase is constant over each grain. Further, is the vector potential, which we here separate into two parts


where is constant in time and corresponds to a uniform magnetic field perpendicular to the array, and is spatially uniform, but time dependent and is necessary to include in order to describe temporal fluctuations in the average electric field , when periodic boundary conditions are used Kim (); Marconi (). Local fluctuations in the magnetic field and hence are ignored.

ii.1 Langevin dynamics

The Langevin dynamics represents a phase-only description of the time-dependent Ginzburg-Landau dynamics (TDGL) in uniform magnetic field. The phenomenological equations of motion for and are of local relaxation type


where the time constant is dimensionless and , is the XY model Hamiltonian, and and are Gaussian stochastic white noise processes.

Figure 1: Equivalent circuit model for Langevin dynamics. The cross denotes the Josephson junction with critical current .

These equations can be cast in the form of circuit equations for an electric circuit built up using the elements displayed in Fig. 1, where each site is connected via a resistor to ground:


The sum in the first equation is taken over the set of superconducting grains connected to and is equivalent to a lattice divergence. In the second equation denotes a sum over all links in the system. The vector goes from site to site and is a fixed external current density. The Gaussian white noise terms and , which now can be interpreted as Johnson-Nyquist noise in the resistors , have the correlations


ii.2 RCSJ dynamics

In the RCSJ model each Josephson junction with critical current is shunted by both a resistor and a capacitor , see Fig. 2. The RSJ model is obtained as a special case when setting . We write the total current from site to site as a sum of the super-, resistive, capacitative, and noise currents


where the voltage between grain and is


and the Johnson-Nyquist noise in the resistors has zero mean and covariance

Figure 2: Equivalent circuit model for RCSJ dynamics.

The equations of motion for the phases and the twists are obtained from local current conservation at each grain


Here we introduced a “transport current” , which is a static, divergence free current distribution, satisfying any external boundary conditions, but is otherwise arbitrary. One may, for instance, connect some nodes to fixed external current sources or sinks. In the present work we will usually use periodic boundary conditions instead, together with the condition of a fixed average current density in the system,


For model purposes we may define a local current density on the links of the lattice as


which directly leads to Eq. (20) when averaged over the system.

ii.3 Lattice structure

We are interested in modeling both regular and random granular superconductors. At low magnetic fields the vortex separation is large compared to the granules and the lattice structure does not matter much. In this regime the models approximate bulk superconductors. In the opposite limit of high magnetic fields the lattice structure is important as frustration effects strongly influence the transport properties. Note that the formulation of the models defined above is independent of lattice structure. We will limit ourselves to two dimensions in the present work. In our simulations presented below in Sec. VI we consider not only square lattices, but also random lattices appropriate as models of random granular superconductors. These irregular lattices are constructed by generating a set of randomly distributed points with unit density on a square and connecting nearest neighbors by Delaunay triangulation. To control the randomness we use the parameter , which is the shortest allowed distance between any two points in the system. Large values of will thus create a relatively ordered lattice structure, whereas lattices with small values are more disordered. For example, a given value of corresponds to a distribution of grain size diameters with standard deviation , respectively. Examples can be seen in Fig. 15. For the random lattices we use two different models, one where the critical current of every junction is set to a constant and one where the critical current is proportional to the contact length between the grains, (in the RCSJ case we also take and ). This length is simply the distance between the points in the dual Voronoi lattice, see Fig. 3.

Figure 3: (Color online) A part of a random lattice with the length and the vector defined. The green lines make up the direct Delaunay lattice and the orange structure is the dual Voronoi lattice.

ii.4 Transport coefficients

In the models defined above the only nonzero transport coefficients are the Nernst signal , the diagonal thermal conductivity and the electrical resistivity . These may be obtained from equilibrium correlation functions using Kubo formulas Kubobook (); Luttinger1964 (). While relates to a mechanical perturbation, and give the response to a nonmechanical property, namely a temperature gradient, and the applicability of a standard Kubo formula is not evident. Nevertheless, we show in the next section that the transport coefficients can be expressed in standard form as


where is the area of the system. In these equations is the average electric field in the -direction and


is the average heat current in the -direction, with and . This form of the heat current, which is one of our main results, will be explicitly derived in the next section for the specific models under consideration.

Iii Heat current

In order to calculate thermoelectric effects and the heat conductivity an expression for the heat current is needed. Several microscopic derivations of the heat current in superconductors has been given in the literature Schmid (); Maki1968 (); CaroliMaki (); Hu (). The presence of a magnetic field yields a complication in that the total energy and charge currents consist of magnetization currents in addition to the transport currents Obraztsov (); CooperHalperinRuzin ().

Rather than relying on these microscopic expressions we derive here, within the framework of the models defined in Sec. II, an expression for the heat current that can be used consistently in our calculations. Below we will first derive the energy current by considering the continuity equation for the energy density. This gives the heat current after subtracting the nondissipative energy transport, except for a magnetization contribution, which is calculated separately in Sec. III.1.3.

Since we are interested in the response of the system to an applied temperature gradient, which is a nonmechanical thermodynamic force, the standard derivation of a Kubo formula does not hold. One approach is to follow Luttinger Luttinger1964 (), and introduce a “gravitational” field, which couples to the energy density, and then proceed with the calculation of the response to perturbations in that field. For the present models, however, there is an alternative route. With the dynamics governed by stochastic equations, Eqs. (12-13) or (19-20), in which the temperature enters via the strength of the stochastic noise, it is possible to introduce local temperature variations, such as a temperature gradient, and calculate the resulting response. This calculation, done in Sec. III.2 gives us both the Kubo formula Eq. (22) and the heat current Eq. (25). Note that in this setting, each point in the model, or more precisely, each resistor in the circuit, is in contact with a local heat bath. For finite gradients one cannot expect the heat current to be automatically conserved, since the resistors act as sinks and sources. For an infinitesimal thermal gradient the heat current will, however, be conserved on average. Alternatively, one could adjust the local temperatures to make sure that heat is conserved on average also for finite temperature gradients. Such self-consistent temperatures have been employed in studies of heat conductivity in harmonic lattices Bonetto2004 (). For the problem at hand, however, the temperature profile is determined by the total heat transport, including phonons, etc, so an externally imposed temperature gradient is probably more realistic. For the linear response the form of the profile should not matter as long as it is smooth. We will use a linear temperature profile below.

iii.1 Heat current from continuity equations

In this section we will derive the heat current expressions for granular superconductors described by Langevin and RCSJ dynamics, as used in our simulations. First, however, it is useful to discuss briefly the continuum formulation. Starting from the thermodynamic relation


for a superconductor in a magnetic field and electric field , where and are the entropy and energy densities, the chemical potential, and the density of charge carriers (with charge ), one obtains for the heat current density


where is the electric transport current density. The total energy current density is the sum of two parts, a nonmagnetic part and Poynting’s vector,


and the transport heat current density therefore also has two contributions LarkinVarlamov (); Ussishkin (); MukerjeeHuse (); Hu ()


where is the magnetization. The latter part can be rewritten as


where is the magnetization current, and the electric potential. The first term on the second line is purely rotational and will not contribute to the heat transport when integrated over the system. The second term combined with the nonelectromagnetic part is the standard expression , where is the electrochemical potential. The last piece is with our gauge choice [Eq. (8)] spatially constant. Therefore it cannot be uniquely determined via the continuity equations below, but will have to be added separately later.

Note the dual role played by the subtraction in Eq. (27). It contains the subtraction (but in a gauge invariant way). At the same time it subtracts the magnetic field energy transported by the vortex current, since can be viewed as a magnetic contribution to the vortex chemical potential, while the electric field is proportional (but transverse) to the vortex current . In principle one should also subtract the nonelectromagnetic energy transported by vortices . In the present models, however, . In fact, the chemical potential for the charge carriers is also zero. We will now derive the analogous expressions for the discrete models.

iii.1.1 Langevin dynamics

It is instructive to study a slightly more general model which includes the charging energy of the grains, described by a circuit model with capacitors added in parallel with the resistors to ground. This modification, besides being more general, makes the derivation more physically transparent, while leaving the final results unchanged. The total energy for this model can be written


with being the gauge invariant phase difference between site and and the voltage at site to ground. With the site we associate a local energy


The time derivative of this is


In the first term we may identify . The last term contains the current through the capacitor , which is eliminated using Kirchhoff’s law


where is the current through the resistor to ground, including the noise current. This results in


This is on the form of a continuity equation for the local energy, since the sum over neighbors is the lattice analogue of a divergence, and the source term represents the dissipated work done by the system on the environment. The term involving the vector potential will be cancelled when the magnetization contribution is dealt with in Sec. III.1.3. The energy current is identified as


and by subtracting the transport current we obtain the heat current


for Langevin dynamics (excluding the magnetization contribution).

iii.1.2 RCSJ dynamics

As in the Langevin case it is convenient to add a capacitance to ground to the usual RCSJ model. The energy for such a model is


where is the voltage across the junction between site and . This implies a local energy of the form


with a time derivative


The last term is, as in the Langevin case, eliminated using current conservation


where the supercurrent , the current through the shunting resistor (including the noise) is , the parallel capacitance current , and the current through the capacitance to ground . We get


and by introducing the total current flowing from to , and rearranging


These terms have interpretations similar to the Langevin case. The second term on the left hand side is the lattice divergence of the energy current, the first term on the right hand side will be cancelled by the magnetization contribution, and the last one represents the work done on the environment. The energy current is thus


and the heat current


for RCSJ dynamics (again excluding the magnetization contribution). The result is very similar to the Langevin case, except that the total current appears instead of just the supercurrent.

iii.1.3 Magnetization contribution to the heat current

The models defined in Sec. II are formulated in the limit where fluctuations in the vector potential are completely suppressed, except for the uniform part . Even so, the latter will, perhaps surprisingly, contribute to the heat current. We will derive this contribution in a more general setting where local fluctuations in the vector potential are allowed. To write down the magnetic energy, we first split the total current flowing between sites and into transverse and longitudinal parts


(this is the lattice analogue of writing a vector field as the sum of a gradient and a curl). The variables are defined on the plaquettes of the lattice, i.e., on the dual lattice, and are often referred to as loop currents, see Fig. 4. The remaining part is the loop current of the loop , which can be nonzero in the presence of resistors and/or capacitors to ground. Without loss of generality we set , so that the heat and energy currents coincide. The loop currents can be used to express the magnetic fluxes through the corresponding loops, via the self- and mutual inductances of the equivalent circuit diagrams (Fig. 12),


where we have chosen our gauge such that , where denotes the ground, and is the applied flux. The sum over is a sum over adjacent plaquettes separated by the link , with and defined by Fig. 4. With this the magnetization energy is


To simplify the argument we only included the self inductances of the loops connecting the ground in the last term.

We may associate a local energy with the plaquette , whose time derivative is


The first term in the last line cancels exactly the corresponding contribution in Eqs. (III.1.1) and (43). The summation in the last term represents the energy flowing into the plaquette from the adjacent plaquettes, and we can therefore identify the magnetization contribution to the energy current and hence the heat current


It is not hard to see that this result holds also in presence of an electric transport current, if the :s are defined as in Eq. (46). Equation (51) is the discrete analogue of the last term of Eq. (III.1).

Figure 4: (Color online) On each plaquette and , defined with respect to the direct lattice points and , there are loop currents and . The difference between these loop currents is the transverse part of the total current on the link from to .

iii.1.4 Full heat current

Averaging Eq. (37) or (45) over the system and adding the magnetization contribution Eq. (51) gives the full average heat current density in the -direction


where is the -component of the difference vector from site to , , and the plaquette centers (the vertices of the Voronoi graph). In the limit where is spatially uniform the last term simplifies to , where is the average magnetization density [ is the area of the (Voronoi) cell ], in agreement with Eq. (III.1). A more practically useful form is obtained below in Sec. III.2.4 by expressing the loop currents in terms of the total currents , which results in the expression Eq. (25).

iii.2 Heat current from Kubo formula

It is useful to see how the heat current enters the linear response arising from an applied temperature gradient via a Kubo formula. We will do this starting from the stochastic equations of motion of either Langevin [Eq. (12)] or RCSJ [Eq. (19)] dynamics. The temperature enters these equations only via the strength of the Gaussian noise correlation function, which now gets a spatial dependence 111In reality the parameters of the model, e.g., , are going to be temperature dependent and therefore also spatially dependent in a temperature gradient. However, this dependence only changes the equilibrium distribution of the model and does not give rise to a thermodynamic force. We can therefore ignore this temperature dependence in the derivation of the heat current. is considered a weak perturbation, and we are interested in calculating the response of some dynamical variable to such a perturbation. We assume open boundary conditions in the -direction, i.e., in the direction of the applied temperature gradient, while in the transverse -direction they can be arbitrary. With no net electric current flowing through the sample the heat current is equal to the energy current.

For concreteness we consider the response of the electric field perpendicular to the temperature gradient and the magnetic field, which gives the Nernst signal . Our goal is to express this linear response using a Kubo formula


where then can be identified with the heat current (the overbar denotes a spatial average), and is the system size.

To this end is convenient to reformulate the problem using functional integrals Janssen1976 () and write ensemble averages as


where is the statistical weight of a given realization of the stochastic process . The Jacobian in Eq. (54) comes from the transformation from the Gaussian noise to Janssen1976 (), but since it plays no role in the following it will be dropped henceforth. The linear response of to a temperature gradient can now be expressed (assuming time translation invariance) as




We are interested in the response to a static perturbation in a stationary state given by


Despite the similarity with Eq. (53), we cannot immediately identify with at this stage, because of their different symmetry properties, an issue we now digress into.

iii.2.1 Symmetry relations of transport coefficients and correlation functions

Before continuing the derivation of the heat current from the Kubo formula we discuss some general properties of linear response in classical statistical systems. It is well known that time reversal symmetry implies symmetry relations among the transport coefficients and among correlation functions. In the presence of a magnetic field the time reversal operation should reverse also that. Consider now the correlation function entering Eq. (55), which in general also depends on the magnetic field . The electric field is even under time reversal, while is neither even nor odd. Instead we may split into even and odd parts. By parity symmetry the correlation functions are odd in , so that


This leads to the following symmetry of the response


For causality implies that , hence , so that the response can be written solely in terms of the odd part of as


where is the Heaviside step function. For the Nernst signal we then get


Thus, in order to evaluate the response it is enough to consider the odd part of .

Let us also make the following observation: At finite times the response will contain transients, which we will not be interested in and which do not contribute to the stationary response. Indeed, any contribution to of the form of a total time derivative will contribute a term


to , which vanishes provided is stationary in equilibrium. We will now treat the Langevin and RCSJ case separately.

iii.2.2 Langevin dynamics

For the case of Langevin dynamics the dynamical action obtains from substituting using Eq. (12) in the Gaussian noise probability distribution


resulting in




and correspondingly


Since is even and is odd under time reversal 222 changes sign under time reversal, the odd part of is


Putting we may rewrite this as


The first term is a total time derivative of times a local energy


and will therefore not contribute to static response functions. The remaining part can be identified with the heat current in the Langevin model,


This form of the heat current is directly formulated using the currents and potentials, and therefore simpler to use than Eq. (III.1.4). To show the equivalence of Eqs. (71) and (III.1.4) it is necessary to go through some further steps. Before doing that, however, we consider the RCSJ case.

iii.2.3 RCSJ dynamics

It is convenient to reformulate the equations of motions for RCSJ dynamics Eq. (19) as


where is the longitudinal part of the total current flowing through the links of the lattice, and is the transverse part, with the defined on the dual lattice sites, i.e., on the plaquettes, adjacent to the bond as in Fig. 4. The :s are often referred to as loop currents. Neither of these contribute to the transport current . As in the Langevin case we set , whereby the heat current equals the energy current.

The dynamical action corresponding to Eqs. (7273) can be expressed as


The first term represents the Gaussian distribution of the white noise current , and is a Lagrange multiplier to enforce the constraints Eq. (72). In a temperature gradient the local temperatures are position dependent, .

The functional integration is over the variables , , , , and . Integrating over and we get