text:
# DESY 03-073June 2003 Lattice QCD Calculations on Commodity Clusters at DESY

## I Introduction

Lattice field theory has established itself as an integral part of high energy physics by providing important, non-perturbatively obtained results for many physical observables. It complements standard approaches of theoretical particle physics such as perturbation theory and phenomenology, becoming an indispensable method to allow for first principles interpretation of experimentally obtained data. The aim of lattice field theory is to understand the structure of field theories, to test these theories against experiment and in this way to find physics behind the standard model of which we know that it has to be there, but not at what energy scale it should appear.

Another important aspect of lattice field theory is its high computational needs. A distinctive feature of the numerical computations in lattice field theory is the required performance of several Tflops and this performance is needed “in one piece”. Such a requirement can only be fulfilled with massively parallel architectures having many processors that are connected via a very fast interconnecting network and work simultaneously on the same problem as a single machine. This distinguishes the computational needs in lattice field theory from the farming concepts usually employed in grid computing of experimental high energy physics.

The combination of high computational needs and the motivation for this in high energy physics as described above, makes a conference such as CHEP an ideal place to present the status and the perspectives of lattice field theory. What strengthens the connection even more is that, at least, a data grid is also needed in lattice field theory to exchange expensive data, the so-called configurations and propagators that are generated. An international lattice Data Grid initiative has been started ildg .

The computational requirements with high performance in the multi-Tflops regime and fine grained communication can be realized on commercial supercomputers like Hitachi hitachi  or IBM ibm , with specialized machines such as APE (Array Processor Experiment) ape  or QCDOC (QCD on Chip) qcdoc , or with PC clusters as they are discussed in this article. Although commercial supercomputers can conveniently be used, since they are maintained by the corresponding computer centers, their price is normally very high, the efficiency of lattice field theory code is often not optimal and they have many users from different application fields leading consequently to a situation that a single user will not get much computer time.

Specialized, custom-made machines like APE and QCDOC are cost-effective and offer the best price/performance value for multi-Tflops installations. On the other hand, a lot of work has to be spent by the physicists themselves in order to develop, build and maintain these machines. A compromise between supercomputers and custom-made machines might therefore be commercial PC clusters and in this article111This paper covers three talks jansen2003 ; gellrich2003 ; wegner2003  given during the parallel session on Lattice Gauge Computing at Computing in High Energy Physics, UCSD, La Jolla, USA, March 24-28, 2003. we will concentrate on the experiences we have gained with these kind of machines at DESY. We also refer to evalgroup  for a thorough discussion on the above point.

One observation in lattice field theory is that people gather in larger and larger collaborations. This results first of all from the huge computer resources required. It is aimed to use these resources wisely, not duplicating results and find the best strategies for solving the physics problems. One example of such an effort is SciDAC scidac  in the US. Another one is the Lattice Forum (LATFOR) initiative in Germany latfor . Associated to LATFOR are also researchers in Austria and Switzerland. In this initiative groups from many universities and research institutes who work on lattice field theory combine their efforts to reach physics milestones in different areas of lattice field theory. They try to coordinate their physics program, to develop and share software and share configurations and propagators, which play the role of very expensive raw data for lattice field theory computations.
The research areas of LATFOR are broad and cover

* ab initio calculations of QCD with dynamical quarks

+ Hadron spectrum and structure functions
    + fundamental parameters of QCD, i.e. the running strong coupling $`\alpha _s(\mu )`$ and the quark masses $`\overline{\mathrm{m}}(\mu )`$
    + B-physics
* matter under extreme conditions

+ QCD thermodynamics
    + QCD at non-vanishing baryon density
* Non-QCD physics

+ Electroweak standard model
    + Supersymmetry
* Conceptual developments

+ exact chiral symmetry on the lattice
    + acceleration of the continuum limit
    + non-perturbative renormalization
    + finite size effects
    + algorithm development

It would be too demanding (and it is not the purpose of this article) to discuss these topics in detail. The main target of lattice calculations is certainly QCD and we will therefore give in the next section a –presumably too short– introduction to QCD and discuss a few examples of the results that can be obtained.

## II Physics Motivation

### II.1 Quantum Chromodynamics on the lattice

In lattice field theory books , the continuum space-time is replaced by a $`4`$-dimensional, euclidean grid with a lattice spacing $`a`$, measured in fm. The advantage of this procedure is that now the theory can be simulated on a computer. In order to obtain back the desired results in the continuum, the values of observables obtained at non-vanishing values of the lattice spacing have to be extrapolated in a continuum limit to zero lattice spacing. At this stage, the comparison to experiments becomes possible and a test of the validity of the model considered can be performed.

On the pure computational side we are dealing with two kind of fields. The first one represents the quarks and are given by complex vectors

$$\mathrm{\Psi }(x)_{\alpha ,a,n_f},\{\begin{array}{cc}\alpha =1,2,3\hfill & \mathrm{color}\mathrm{index}\hfill \\ a=1,2,3,4\hfill & \mathrm{Dirac}\mathrm{index}\hfill \\ n_f=1,\mathrm{},6\hfill & \mathrm{flavor}\mathrm{index}\hfill \end{array}.$$
(1)
These fields live on the $`4`$-dimensional space time point $`x`$ which, on the computer, is represented by integer numbers, $`x=(x_1,x_2,x_3,x_4)=a(i,j,k,l),1i,j,k,lN`$. A second set of fields represents the gluons of QCD and are given by SU(3) matrices $`U`$, $`3\times 3`$ complex matrix with unit norm. The fields $`U(x,\mu )_{\alpha ,\beta }`$ carry again color indices through which they interact with the quark fields. The gluon fields live on the links of the lattice that connect points $`x`$ and $`x+\mu `$ in direction $`\mu =1,2,3,4`$. The interaction is described by the action222Of course, in QCD the quark fields are represented as Grassmann variables. We discuss here the bosonized form of the action as it is used in simulations.

$$S=\overline{\mathrm{\Psi }}M^1\mathrm{\Psi }.$$
(2)
The action of eq. (2) requires the inverse of the so-called fermion matrix $`M`$, or, to be more specific, the vector $`X=M^1\mathrm{\Psi }`$. Without giving the exact definition of the matrix $`M`$, the problem is that the fermion matrix is high dimensional $`O(10^6)O(10^6)`$ and the numerical solution of the linear set of equations

$$MX=\mathrm{\Psi }$$
(3)
employing such a matrix is clearly very demanding. It helps, however, that the matrix is sparse. For such a case a vast literature for solving eq. (3) exists saad . What is important is that the algorithms that can be employed are self-correcting and many of them can be proven to converge in at most a number of steps that corresponds to the dimension of the matrix. Thus we are left with a well posed and regular numerical problem. Of course, in practice the number of iterations is much smaller and typical numbers of iterations to solve eq. (3) are $`1001000`$. Note that in each of these iterations the matrix $`M`$ has to be applied to a vector of size $`O(10^6)`$.

In order to give a feeling about the computational demand, let us give an example. Let us consider a lattice of size $`32^364`$ as is realistic in todays calculations in the quenched theory, where internal quark loops are neglected. Then we would need to solve eq. (3) twelve times per configuration for each color and Dirac component. Each solution needs $`O(200)`$ iterations and we want to perform this on typically $`O(1000)`$ configurations. Since one application of the matrix $`M`$ on such a lattice needs $`2.8`$ Gflop, we are left with approximately $`6.6`$ Pflop for obtaining only one physical result for a single set of parameters, i.e. the bare coupling and the bare quark masses. On your standard PC which might run with $`500`$ Mflops sustained for this problem, you would hence need about five months. In order to reach control over the finite size effects, the chiral extrapolation and the continuum limit, simulations at many values of the bare parameters have to be performed.

The numbers above hold for the quenched case, where the quark fields are left out as dynamical degrees of freedom in the simulation. If they are included, the cost of the simulations becomes at least a factor $`100`$ more and a single physical result would need about $`40`$ years on your PC. Clearly, better computers are needed such as the ones discussed in the introduction.

Of course, relying only on progress in the development of computers would be too risky and not enough. Lattice field theory has seen a number of conceptual improvements oaimprove ; chiral  in the last years that allowed to accelerate the simulations themselves. In addition, many improvements in the algorithms used were found. Although each algorithmic improvement by itself was only a relatively small step algoreview , all in all a factor of $`20`$ acceleration through algorithm improvement alone could be obtained in the last $`15`$ years. Still, the development of machines were much faster in this period.

The situation is illustrated in fig. 1 (taken from a LATFOR paper). Here we show the speedup obtained relative to the status in the year 1987. This year is special in that at that time the first exact algorithm for simulations of dynamical fermions was developed and used hmc . We see in the figure that since then the speedup resulting from better computer techniques, despite the impressive improvements by algorithmic developments, is orders of magnitude larger than the algorithm improvement. The figure also shows that special hardware, in this case the APE computer, shows the same scaling law as commercial supercomputers, in this case the CRAY. Other computers and their performance, relative to the status in 1987, used throughout the world for lattice field theory is also shown.

### II.2 Some selected physics results

The machine, algorithmic and conceptual developments in lattice field theory allowed to compute a number of important physical quantities in the last years, despite the aforementioned very large computational requirements. An important aspect of the developments is that we understand not only the statistical errors of the numerical calculations, but also the systematic ones. For a number of quantities fully non-perturbatively renormalized results in the continuum limit could be obtained. The only restriction of these calculations is that they are still done in the quenched approximation. However, there is no particular reason why the calculation should not be doable in a completely analogous way for the full theory. The advent of the next generation of machines in the multi-Tflops regime will then open the door to this exciting perspective.

Let us just discuss two examples of physics results to illustrate what we just said. One is the running coupling rainerhart  and the other are moments of parton distribution functions in deep inelastic scattering. In a quantum field theory such as QCD, there is a steady generation of virtual particles that shield (or anti-shield) the charges of the elementary particles. This leads to renormalization effects. By changing the energy at which experiments are performed, the scale at which we look at the, say, color charge is altered, and, correspondingly the value of the charge itself depends on this energy scale.

This scale dependence (the running) of the coupling can be computed in lattice field, starting from the QCD Lagrangian alone, without any further assumptions. The trick is to use a suitable lattice renormalization scheme, the so-called Schrödinger functional (SF) scheme, that is defined in a finite volume and hence ideally suited for numerical simulations. By going to very high energies, contact to perturbation theory can safely be established and renormalization group invariant quantities can be extracted. In the case of the running strong coupling this corresponds to the $`\mathrm{\Lambda }`$-parameter of QCD. The advantage of the knowledge of renormalization group invariant quantities is that they can be translated to any preferred renormalization scheme. In this way it becomes possible to translate results for the running coupling from lattice simulations to continuum results in the, say, $`\overline{\mathrm{MS}}`$-scheme as it is conventionally used in perturbation theory.

In fig. 2 we show an example of such a calculation. The results are already in the continuum. They cover a broad energy range and are very precise (the error bars of the simulation points are well below the size of the symbols). In the plot, also a comparison to perturbation theory is shown and a good agreement down to surprisingly small energy scales is found.

Another example is a moment of a parton distribution function as they can be extracted from global analyses of experimental data. Such moments can be expressed as expectation values of local operators and are hence computable in lattice simulations. The renormalization procedure of such moments follow the general strategy of using the finite volume SF renormalization scheme discussed above for the running coupling. We show in fig. 3 an example of the continuum limit of the first moment $`x`$ of a twist-2, non-singlet operator in a pion structuref . In the plot two different lattice formulations of QCD were used. It is reassuring that in the limit that the lattice spacing is sent to zero both unphysical lattice versions of QCD extrapolate to the same number.

As a final result for this study it is found that the lattice gives a preliminary value of $`x=0.30(3)`$ while the experimental number reads $`x=0.23(3)`$. This lattice number has to be taken with care since it is obtained in the quenched approximation. But, again, there is no other reason than missing computer power to repeat this calculation also for the full theory. If in this case the results just mentioned were stable, a really interesting situation would have emerged.

The results that we have just discussed were actually obtained on APE machines. There are, however, a number of physics problems, where PC clusters were used extensively. The simulations concern in particular the recently discovered chirally invariant formulations of QCD on the lattice. Examples for results on the PC clusters that are installed at DESY are given in clusterresults . The next sections are devoted to a discussion on PC cluster systems that are installed at DESY.

## III Commodity Clusters at DESY

### III.1 Conceptional considerations

Sufficient computing power to perform Lattice QCD (LQCD) calculations as described above can obviously not be drawn from a single processor. The solution is to parallelize the physical problem in order to concurrently deploy many CPUs. As a consequence massive intercommunication between the CPUs is needed.
Typical High Performance Computing (HPC) or Supercomputing applications are characterized by high demands on:

* CPU (especially Floating Point Unit (FPU)) performance,
* memory throughput,
* interconnectivity (bandwidth and latency) between nodes.

One example for so-called supercomputers are Symmetric Multi-Processing (SMP) machines with up to hundreds of processors. They appeared as single machines and are optimized for memory access and interconnectivity between the CPUs.
In the LQCD area custom-made special purpose machines such as APE (Array Processor Experiment) ape  or QCDOC (QCD on Chip) qcdoc  have been developed.

In the last years PCs, exploiting processors with competitive computing power, hit the commodity market. Concurrently, modern network technologies have achieved performances, which allow for high-speed low latency interconnects between PCs. These developments paved the way to build PC clusters clusters .
Commodity clusters draw computing power from up to hundreds or even thousands of in principle independent PCs with one or two CPUs, called nodes. Those clusters benefit from their scalability and the possibility to deploy components of the commodity market with good price/performance ratios. Interconnectivity is provided by exploiting modern network technologies. Such a cluster must:

* deliver sufficient CPU (FPU) performance and memory throughput,
* provide good connectivity between CPUs (bandwidth as well as latency),
* be scalable,
* be reliable,
* incorporate tools for easy installation and administration,
* provide a usable software environment for the applications,
* be connected to backup and archiving facilities,
* fit boundary conditions such as space, cooling and power supply capacities.

In clusters, a set of main building blocks can be identified:

* The computing nodes which actually provide the computing power, optionally with local disk space,
* a high-speed, low latency network for the parallelized physics application,
* an auxiliary network to remotely control and administer the nodes,
* a host system for login, compiling, linking, batch job submission, and central disk space,
* optionally, a slow control network, e.g. based on a field bus.

A schematic view is shown in fig. 4. The high speed network can either be organized as a switched network (e.g. the DESY clusters using Myrinet-switches) or by a n-dimensional mesh to allow for nearest-neighbor communication, see fig. 5. In Budapestcluster  the network is organized as a $`2`$-dimensional GigaBit-Ethernet mesh.

In scientific computing Unix-like operating systems have always played the dominant role. The development of Linux along with the triumphal procession of PCs into the scientific world made Linux-PCs the systems of choice in most universities, physics institutes, and laboratories. In order to actually operate Linux-PC clusters, further system aspects must be taken into account:

* Installation and administration of the operating system Linux,
* security issues (login, open ports, private network),
* user administration,
* application software installation,
* backup and archiving issues,
* monitoring and alarming.

### III.2 Implementation

The recently improved PC architectures are well suitable for Supercomputing by exploiting:

* Increasing CPU clock rates following Moore’s Law now extending the $`2`$ GHz border,
* larger caches at full processor speed,
* vector units (SSE, SSE2),
* cache pre-fetch,
* fast memory interfaces,
* PCI-Bus at $`66`$ MHz/$`64`$-bit,
* high external bandwidth (Myrinet, GigaBit-Ethernet), partly with very low latency (Myrinet).

#### III.2.1 Hardware

At both DESY sites in Hamburg and Zeuthen commodity clusters are operated since January 2002 and December 2001 respectively. See tab. 1 for the set-ups.

The cluster nodes as well as the host system are equipped with high-end commodity components (see tab. 2).

Intel Pentium 4 Xeon processors xeon , which became available at the end of 2001, showed much enhanced performance compared to Pentium III Tualatine CPUs due to features such as vector units, and larger caches (see the section IV on benchmarking results). The SuperMicro motherboard P4DC6 with the Intel i860 chip-set was the only possible combination until the end of 2002. It provides PCI-Bus support. For the memory Rambus modules were chosen. The communication between the nodes in the physics application is done by means of Myrinet myrinet . It provides bandwidths up to $`240`$ MB/s with very low latencies in the order of a few $`\mu `$s. For administration purposes and to actually submit jobs and copy data the nodes are interconnected via Fast-Ethernet in a private subnet. The host system is connected to the central switch by a copper GigaBit-Ethernet link and has a separate fiber GigaBit-Ethernet link to the outside. Each node consists of a PC with two CPUs and is housed by a 4U rack-mounted chassis. Up to $`9`$ nodes ($`18`$ CPUs) are installed in a cabinet.

#### III.2.2 Software

The basic installation plus most of the software support was purchased with the hardware from the German company MEGWARE megware .
For all software related aspects the host system is used as a server for the nodes.
For the Linux operating system a S.u.S.E. 7.2 distribution suse  was chosen, which contains the kernel version 2.4.17 with SMP capabilities. Temperature and fan sensors are read out by the kernel module lmsensors lmsensors .
The nodes are booted via DHCP, TFTP, and Intel’s intel  Pre-boot Execution Environment (PXE).
The nodes are operated in a private network (192.168.1.0) behind the server. For security reasons, external user login is only possible to the host system via ssh. Individual login from the host system to the nodes can be done over the auxiliary network by means of rsh. Users are registered on the host system which exports the home directories (/home) and a data partition (/data) to the nodes. For the user administration standard Unix tools are used (useradd). The necessary files (/etc/groups, /etc/passwd, /etc/shadow) are manually distributed to the nodes.
Source code is compiled and linked on the host system. Software is mostly written in C/C++ but also compilers for Fortran77/90 are required. On the host system in addition to GNU compilers of the Portland Group pgi , KAI kai , and Intel intel  are available.
The parallelization of the computation is done within the application by means of the Message Passing Interface (MPI). Since MPI is running over Myrinet, a special library which uses low level Myricom communications is installed (MPICH-GM mpich-gm ). For the Myrinet network a static routing table is used.
Zeuthen uses the open source Portable Batch System (OpenPBS) for job submission. In Hamburg nodes are manually distributed to users on good-will basis.
Time synchronization is done via XNTP. The nodes synchronize with respect to the host system which gets the correct time from DESY’s central server.

Backup and archiving are basically different items: Regular backups are done to provide security against loss of system data and home directories. Under normal conditions backups will never be retrieved. For the Hamburg cluster DESY’s standard backup environment based on IBM’s Tivoli Storage Manager (TSM) is used. It automatically creates incremental backups of the disk of the host system and a regular basis. At DESY Zeuthen a copy of the host system’s disk is stored on a second disk.
Archiving tools allow users to arbitrarily store and retrieve large amounts of data. DESY uses dCache which provides a simplified and unified tool to access the tertiary storage dcache . It provides a unique view into the storage repository, hiding the physical location of the file data, cached or tape only. Dataset staging and disk space management is performed invisibly to the data clients. Currently around $`1`$ TB of data are stored. The archiving system is mainly used to store temporary check-points and final results of long time computing jobs, so-called configurations, which can be used for further analyses.
The backup and archiving scheme is shown in fig. 6.

LQCD calculations require the availability of all nodes dedicated to the problem for the entire run-time of the job. The failure of only one node spoils the entire job. Therefore, stability of the cluster in terms of availability and sustained performance is a crucial credential of the clusters. Usually not all nodes are used for one single job. Many calculations are done on a limited number ($`28`$) nodes. This allows to (manually) restart jobs using the check-points on different nodes in case of failures. In order to use the resources of the clusters most efficiently a well-defined monitoring and alarming scheme is needed.
The company MEGWARE delivered with the software installation a monitoring package called clustware. It provides a snapshot of all relevant properties of the cluster nodes in a graphical user interface, including CPU usage, load, I/O, and temperature. A long term history ($`>1min`$) is not shown.
Alternatively, a DESY in-house development called ClusterMonitor (CluMon) clumon  is in use. Every node runs a simple Perl-written daemon which periodically dumps status information of relevant node properties such as uptime, load, CPU usage, memory usage, swap usage, and all temperatures into a node-specific file on the host system. The host system runs an Apache web server, which allows to remotely access status information of the cluster with any web browser. History of the quantities is kept by means of the MRTG package mrtg  and is available from the web page. In addition, CluMon provides alarming by e-mail based on the time period since the last update of the status file.

### III.3 Operational experiences

The two clusters in Hamburg and Zeuthen are located in the computer centers. The clusters are operated and administrated in cooperation by members of the computer centers and the DESY theory and NIC groups.
Of around twenty registered users per cluster just a handful can be classified as Power Users, running regularly resource consuming jobs (see tab. 3).

The two clusters in Hamburg and Zeuthen deploy in total $`52`$ dual-CPU PCs: $`32`$ nodes (Hamburg), plus $`16`$ nodes (Zeuthen), plus $`2`$ spare nodes, plus $`2`$ servers. During the almost $`17`$ month of operation quite a number of problems occurred, so usually not all nodes have been available at all times. Tab. 4 lists all major hardware problems.

Taking into account that the clusters exploit commodity hardware components and offer considerably better price/performance ratios than big mainframe SMP-machines, failures of certain components such as disks and power supplies were expected. The stability of the PC hardware after replacing the obviously systematically misbehaving IBM disks was reasonable. Some annoyance was caused the repeating failures of the Myrinet interface and line cards, which was also seen at Fermilab holmgren  and traced back to broken optical receivers. Nevertheless, the general opinion of the users on the performance and the stability of the clusters is very positive.

### III.4 Future developments

The considerations in section II require cluster sizes of $`O(1000)`$ nodes to approach the Tflops regime. Even more, in order to actually deliver a few Tflops sustained for hours, days or even weeks, all nodes would need to run at the same time. As discussed earlier the failure of just one node would spoil the entire calculation. Accepting the experiences so far, this seems not be possible within the current concept333$`1`$ broken node per week in a $`50`$ node cluster is equivalent to a maximal lifetime of a complete $`1000`$ node cluster of $`8`$ hours..
Recent tests have shown that GigaBit-Ethernet might be an interesting alternative to Myrinet. Benchmark showed bandwidths of $`2\times 1`$ Gbit/s bidirectional with special switches which would imply a much better price/performance ratio than Myrinet. Since GigaBit-Ethernet is widely used now, one might also expect more stability and reliability compared to the niche product Myrinet (see section IV).
Disks –even SCSI– are the most likely components to break in PCs. Though the replacement of disks is easy, the affected node is down and needs to be re-installed completely afterwards. This could be avoided by running the nodes disk-less. Such a concept would require a stable and reliable server which could be achieved by setting up one or more RAID-systems in a tree-like architecture to distribute load. The server could also provide the boot image. In another scenario booting could be done from a local block device such as an EPROM or a memory stick.
The current set-up relies on one single server which exports /home and /data directories to the nodes. It also serves as a login host for the users and is used for code development, compilation, linking, and job submission. This machine is clearly a single point of failure. In a bigger system one would opt for redundancy in the server arrangement by distributing different functionalities to different machines. In particular a separate file server for the exported directories is needed.
Space, power consumption, and cooling will become a major issue when planning for thousands of nodes. Recent developments of so-called blades place the motherboards vertically to improve the air-flow for cooling in order to increase CPU densities.
Software installation, administration, and monitoring of thousands of nodes is a challenge which requires a very careful choice of appropriate tools. Remote administration could be enabled by exploiting the serial consoles of the PCs. They could be connected to a dedicated terminal server or –in a much cheaper scenario– subsequently from node to node.

## IV Benchmarks

The hardware of commodity PCs has been extremely improved over the last years. The performance increase inside the CPU is due to higher clock rates and enhanced building blocks (e.g. SSE1/SSE2 instructions) following Moore’s Law which predicts a performance doubling all $`18`$ month. Moore’s Law gives a technology estimation of mainly the CMOS density or number of transistors which can be integrated on a chip of a given size. It does not work well for the other interacting PC components like the memory interface and external busses. On the other hand also a big step forward in the development of fast memory architectures like Rambus and DDR RAM and a series of high bandwidth PCI bus based interconnects like Myrinet and QSNet is going on. Therefore PC clusters are becoming more and more attractive for classical high performance parallel computing and therefore also as a hardware basis for LQCD applications.

### IV.1 Benchmark systems

Apart from the DESY clusters described above, the following systems were used in order to test the ability of PC clusters for LQCD applications:
Mellanox: Blade dual Pentium 4 Xeon cluster connected via Infiniband, running MPICH for VIA/Infiniband with patch from Ohio State University mellanox ,
ParTec: Dual Pentium 4 Xeon cluster connected via Myrinet running ParaStation MPI partec ,
MEGWARE: Dual Pentium 4 Xeon cluster connected via Myrinet running MPICH-GM from Myricom mpich-gm ,
Leibniz-Rechenzentrum Munich: (single CPU tests) Pentium 4 and dual Xeon PCs with CPUs with clock rates between $`2.4`$ and $`3.06`$ GHz,
University of Erlangen: GigaBit-Ethernet dual Pentium 4 Xeon cluster.

### IV.2 Benchmarks and results

Representative benchmarks for the evaluation of different PC systems have been developed. Already in the year 2000 a first benchmark of M. Lüscher (CERN) has shown the potential in using the SSE1 and SSE2 instructions for the Wilson-Dirac operator Luscher2001 . This program takes heavily advantage of the Pentium 4 memory-to-cache pre-fetch capabilities and the SSE registers and instructions which are implemented by using assembly in-line code, compatible to the gcc and Intel compilers. Fig. 7 shows on the left hand side the performance gain of the highly optimized $`32`$-bit and $`64`$-bit Dirac-Operator kernel which linearly follows the evolution of the CPU performance expressed by their clock rate. The value of $`1.5`$ Gflops for the $`32`$-bit implementation or $`0.8`$ Gflops for the $`64`$-bit implementation respectively was unexpected high for a PC in the year 2000 and encouraged groups working on LQCD algorithms on PC clusters also to use the Pentium 4 capabilities to improve their algorithms on PC clusters.
The Wilson-Dirac operator Benchmarks are accompanied by two tests called add\_assign\_field (similar to the BLAS daxpy) and square\_norm which are representing the linear algebra part of the benchmark. Both parts are strongly memory bound which means that they cannot benefit from the SSE-environment. This results in a relative small improvement shown on the right hand site of fig. 7 which also gives an impression of the slowly evolving memory interface architectures since the introduction of the dual channel $`800`$ MHz Rambus.
Another version of such a single node benchmark was developed by M. Hasenbusch (DESY) for the even-odd preconditioned Wilson-Dirac operator cluster2002 . Meanwhile (using recent FSB800 based PCs equipped with a $`3.06`$ MHz Pentium 4) a performance of about $`2.6`$ Gflops for the $`32`$-bit implementation and about $`1.4`$ Gflops for the $`64`$-bit implementation can be observed.
To evaluate the behavior of PC cluster interconnects a $`1`$-dimensional parallel even-odd preconditioned Dirac Operator Benchmark on a $`2\times 16^3`$ lattice (also written by M. Hasenbusch) was used (see Appendix). The aim of the parallel benchmark was to compare different parallel PC based architectures against each other rather than achieving the best performance for a given system. Fig. 8 shows the results on clusters with different numbers of nodes. In addition to the CPU power and the memory interface the throughput of the external PCI-Bus depending on the given chip-set and the interconnecting interface itself are dominating the entire performance. Both early Intel Pentium 4 Xeon based clusters at DESY are using the i860 chip-set which came with a relative poor $`33`$ MHz/$`64`$-bit PCI-Bus performance. A bus-read (send) of $`227`$ MB/s and a bus-write (recv) of $`315`$ MB/s of maximal $`528`$ MB/s and expected $`450460`$ MB/s was measured. This ends up in an external unidirectional bandwidth of about $`160`$ MB/s of maximal $`240`$ MB/s. In the bidirectional case we measured $`(90+90)`$ MB/s. Meanwhile more advanced chip-sets like the E7500 were available for the Pentium 4 Xeon which are providing a PCI throughput close to the expected numbers. The influence of the chip-set is dominating the results shown in fig. 8, whereas in the left hand side on CPU per node and an the right hand side two CPUs per node communicating via shared memory MPI are used.
Infiniband is a new promising communication technology especially designed for cluster interconnections fig. 9. Beside the high throughput shown in fig. 10 the latency at relative small buffer size in the order of $`2`$ kB is significantly higher then using Myrinet which could be an advantage for applications using small local lattices. In fig. 11 the performance of a $`4`$ node partition of the DESY Myrinet Pentium 4 Xeon cluster is compared to the performance of a corresponding Infiniband cluster from Mellanox using a $`2`$-dim parallel Wilson-Dirac Operator Benchmark developed at DESY Hamburg by M. Lüscher. Due to some problems using the assembly in-lines within the 2.96 gcc compiler coming with the RedHat Linux system on the Mellanox cluster during the short time in which the cluster was available for testing a code without the SSE optimizations was used. The Infiniband cluster performed approximately $`1.8`$ times better in the $`32`$-bit case and approximately $`1.6`$ times better in the $`64`$-bit case as the Myrinet cluster.
The $`1`$-dimensional (non-optimized) parallel even-odd preconditioned Dirac Operator Benchmark was modified to test whether an improvement using non-blocking MPI communication functions can be achieved. No effect was seen in using MPICH over GM provided by Myricom. Tests on a $`4`$ node cluster which runs the MPI version of the ParaStation software results in an performance increase of about $`20`$% (see tab. 5).

Fig. 12 gives a summary of the communication behavior of the different architectures resulting from the parallel benchmarks. The efficiency number is the ratio between the performance of the same benchmark with and without communication. Included is also a test on a $`4`$ node GigaBit-Ethernet cluster connected via a non-blocking GigaBit-Ethernet switch. The efficiency number for GigaBit-Ethernet in the rightmost column of fig. is 12 even better than one can achieve using Myrinet on a system with a chip-set which provides a slow PCI-bus throughput. The positive influence of non-blocking communication support (ParaStation) and fast communication support (Infiniband) is shown by the efficiency numbers.
Compared to the single node benchmarks the results of the parallel benchmarks imply that that the capacity of bandwidth is crucial for the efficient use of PC clusters as a scalable platform for LQCD applications.

## V Conclusions

We have discussed the usage of commodity clusters based on PCs for LQCD calculations. Such installations would not have been considered competitive a few years ago. However, our experience with such kind of machines at NIC and DESY adds further evidence that for problems in QCD that require below, say, 1 Tflops computer power, PC clusters are a valuable and cots-effective tool for computing physics results in LQCD. In Hamburg and Zeuthen clusters with $`64`$ and $`32`$ CPUs are successfully in operation for more than one year.
Investigations using representative benchmarks on the DESY clusters and also other architectures were carried out with promising results.
Applying the SSE/SSE2 (SIMD+pre-fetch) instructions on Pentium 4 like CPUs, the single node performance of the Wilson-Dirac operator is increasing according to the clock rate improvements of those CPUs used in commodity PCs.
The performance of memory bounded parts of the LQCD algorithms, especially the linear algebra routines, depends strongly on the throughput of the memory interfaces. Those interfaces did not show the same level of enhancements as it has been observed and will be expected further in the case of the CPUs.
The performance of the $`1`$-dimensional and $`2`$-dimensional parallel implementations of the Dirac-Wilson operator depends on the behavior of the external interconnects, i.e. is mainly dependent on the PCI-bus throughput coming given by the chip-set and the interface card itself. Results coming from PC clusters consisting of different components have shown an enhancement in both the quality of the chip-sets (e.g. E7500)and the throughput of the communication interfaces (e.g. Infiniband).
Non-blocking MPI communication can improve the performance by using adequate MPI implementations (e.g. ParaStation).

In summary, it might be envisaged, as done by e.g. LATFOR, that heterogeneous computer landscapes will be available to the user with centers that host machines in the multi-Tflops regime, still enabled by specialized machines, and many smaller installations at universities as well as research centers in the few hundred to 1 Tflops range realized by PC clusters.

## Appendix: Discussion of the even-odd benchmark

The benchmark program applies the even-odd preconditioned Wilson-Dirac matrix that is defined on a $`L^3\times T`$ lattice to a spinor-field. The program is implemented in C plus some in-lined SSE2 extensions. For parallelization we have used the MPI message passing library. The code is derived from Martin Lüscher’s benchmark code presented at the lattice conference 2001 Luscher2001 . Even with communication switched off, the present code performs worse then the one of Martin Lüscher ($`579`$ Mflops vs. $`880`$ Mflops on $`1.7`$ GHz Pentium 4). The reason is twofold:

* Less variables that reside in the cache can be reused than in the standard case.
* We have skipped the cache-optimized order of the lattice-points to simplify the parallelization.

### Strategy of the parallelization

The Wilson-Dirac matrix is a sparse matrix. The hopping part of the matrix only connects nearest neighbor sites on the lattice. Therefore, for parallelization, it is natural to divide the lattice in sub-blocks of size $`t\times l_x\times l_y\times l_z`$. Each of the MPI-processes takes one such sub-block. For simplicity we have done the parallelization only in one direction: $`l_x=l_y=l_z=L`$ and $`tn_p=T`$, where $`n_p`$ is the number of processes.

The hopping part of the Wilson-Dirac matrix connects nearest neighbor sites on the lattice. Therefore each application of $`H_{eo}`$ or $`H_{oe}`$ ($`H_{oe}`$ connects even with odd sites and $`H_{eo}`$ vice versa.) the spinor-fields at the right boundary of the left neighbor and the spinor-fields at the left boundary of the right neighbor of each of the processes has to be sent and received.

In the case of even-odd pre-conditioning the spinor-field only resides on the even (or odd) sites. Therefore $`L^3/2`$ spinors have to be sent and received. A single data package has the size

$$24\times 8\times L^3/2\text{Byte}=96\times L^3\text{Byte}$$
(4)
In our blocking implementation, the communication of the data and the computation is performed in a consecutive way. First the spinor-fields are exchanged using the $`\text{MPI}\mathrm{\_}\text{Sendrecv}`$ function. This is followed by the application of $`H_{eo}`$ or $`H_{oe}`$ on the single nodes.

Both times required for communication $`t_{comm}`$ and calculation $`t_{calc}`$ are measured separately. The effective bandwidth is computed as:

$$\text{Bandwidth}=\frac{96\times L^3/2\text{Byte}}{t_{comm}}$$
(5)
The performance per node without communication is computed as:

$$P_0=\frac{t\times L^3/2\times 1392\text{flop}}{t_{calc}}$$
(6)
Correspondingly, the performance including the communication of the data is given by:

$$P=\frac{t\times L^3/2\times 1392\text{flop}}{t_{comm}+t_{calc}}$$
(7)
In the case of the non-blocking communication, we had to divide $`H_{oe}`$ (or $`H_{eo}`$) into a part that only acts on spinors that reside on the local lattice and a part that acts on spinors that reside on the neighbors:

* Initialize send and receive ($`\text{MPI}\mathrm{\_}\text{Isend}`$ and $`\text{MPI}\mathrm{\_}\text{Irecv}`$),
* perform the calculation for the local part of $`H_{oe}`$,
* Wait for the communication to finish ($`\text{MPI}\mathrm{\_}\text{Wait}`$),
* do the rest of $`H_{oe}`$.

###### Acknowledgements.

The authors would like to thank Martin Lüscher (CERN) for the benchmark codes and the fruitful discussions about PCs for LQCD, and Isabel Campos Plasencia (Leibnitz-Rechenzentrum Munich), Gerhard Wellein (Uni Erlangen), Holger Müller (MEGWARE), Norbert Eicker (ParTec), Chris Eddington (Mellanox) for the opportunity to run the benchmarks on their clusters.
The authors also wish to thank the computer centers of DESY Hamburg and Zeuthen.


format: text

text:
# Andreev drag effect in ferromagnetic-normal-superconducting systems

## I Introduction

At a normal-superconducting interface, Andreev reflection causes a conversion of the quasiparticle charge *and* a flipping of its spin—an electron with spin up, incoming towards the superconductor, is reflected as a hole with spin down. This leads to exciting consequences in the transport physics where the focus is put on creating, manipulating and detecting spin currents and magnetization.wol01 ; fie99  When the carriers are spin polarized, the transport properties are strongly sensitive to the relative orientation of the ordered moments of ferromagnetic electrodes,bai88 ; slo89  a feature that has received a lot of interest due to the relatively long spin dephasing times observed in metalsjoh85 ; jed01  and semiconductors.kik97  In hybrid systems containing a superconductor, various aspects of the interplay between Andreev reflection and ferromagnetism have been addressed theoreticallytheoryFM ; bel00  and experimentally.expFM ; expFM2 ; expFM3

The structure of interest here is shown in Fig. 1(a): two ferromagnetic contacts are connected to a normal conductor which is in turn connected to a superconductor. Of particular interest is the case of crossed Andreev reflection between ferromagnetic leads with antiparallel polarization dominating over crossed normal reflection. In such a geometry, with a voltage applied to one of the leads while the other lead and the superconductor are grounded, the injected current from the biased ferromagnetic lead effectively drags along a current from the grounded lead. This effect is easy to understand since the transfer of Cooper pairs into the superconductor requires a pair of electrons with opposite spin. If the injected current consists only of, say, spin-up carriers, the system supplies the spin-down carriers from the oppositely polarized ferromagnetic lead even though no voltage is applied at this contact.

The physics of the Andreev drag effect was already addressed by Jedema and van Wees et al.jed99  for a multi-terminal metallic diffusive conductor. These authors pointed to negative four-terminal resistances as a consequence of the Andreev drag effect. In contrast, much of the recent theoretical discussion has focused on a geometry where two needle-shaped ferromagnetic leads with antiparallel polarization are directly coupled to a superconductor at two spatially separated points.deu00  In this case, even in the case of $`100\%`$ spin polarization, the crossed Andreev reflection probability is strongly suppressed in realistic geometries,cho00  making the experimental detection of the drag effect difficult. In addition, the direct connection of the ferromagnetic electrodes to the superconductor might lead to unwanted modifications of the local properties of the ferromagnetexpFM2  and the superconductor,expFM3  suppressing the drag effect.

In this work, we overcome the above difficulties by inserting a normal multichannelles01  mesoscopic conductor with generic elastic scattering between the ferromagnetic electrodes and the superconductor \[see the sketch of Fig. 1(a)\]. We include arbitrary interface conductances and spin-flip relaxation (see below). The advantage of this scheme is twofold: (i) since the superconductor is only contacted at one point, the drag effect develops to the fullest extent possible, being of the order of the conductances of the contacts between the normal conductor and the ferromagnetic and superconducting electrodes; and (ii) the sources of the competing phenomena—superconductivity and ferromagnetism—are spatially separated by the normal conductor.

We find that the Andreev drag effect is observable when the conductance of the normal-superconducting (NS) interface dominates over the ferromagnetic-normal (FN) contact conductances. The drag effect is reinforced with increasing lead polarization and survives for a considerable amount of spin-flip scattering. Furthermore, we demonstrate that the crossed Andreev reflections have a profound influence on the zero-frequency current–current cross-correlations. We observe that both the thermal noise and the shot noise power measured between the two ferromagnetic leads yield *positive* values.

## II Model

We analyze a hybrid system (shown in Fig. 1) consisting of a normal conductor connected to two ferromagnetic (F1 and F2) and one superconducting (S) reservoirs. For simplicity, the direction of the magnetization of the two ferromagnets is assumed to be collinear and we consider identical FN couplings. Any type of contacts, e.g. diffusive-, ballistic- or tunnel contacts, can be treated. In what follows, we deal with point contacts at the FN (NS) interface with $`N`$ ($`M`$) transversal modes and mode independent transparency $`\mathrm{\Gamma }`$ ($`\mathrm{\Gamma }_S`$). We consider the case with dimensionless conductances (in units of $`2e^2/h`$) much larger than unity: $`g=N\mathrm{\Gamma },g_S=MR1`$, where $`R=\mathrm{\Gamma }_S^2/(2\mathrm{\Gamma }_S)^2`$ is the Andreev reflection probability. Under these conditions, weak localization as well as Coulomb blockade effects can be safely neglected. Moreover, the dwell time of the normal conductor, $`\tau _d`$, is assumed to be shorter than typical inelastic scattering times.

Throughout the paper, we consider the case when the proximity effect inside the normal conductor is suppressed (but the Andreev reflection at the NS interface is taken into account). A strong proximity effect would quench the drag effect as the normal conductor would behave effectively as a superconductor. note  This can be avoided either with a magnetic field (e.g., a stray field from the ferromagnets or an externally applied field) or with a characteristic quasiparticle energy (bias $`eV`$ or temperature $`k_BT`$) much larger than the inverse dwell time, $`eV,k_BT\mathrm{}/\tau _d`$. Nevertheless, to avoid single particle transport, both $`eV`$ and $`k_BT`$ should be smaller than the gap of the superconductor.

The resistance of the system is completely dominated by the resistances of the FN and NS contacts. Moreover, we assume that the physical properties within the conductor are isotropic due to diffusive or chaotic scattering. We can thus consider the normal conductor to be effectively zero dimensional. This restricts our model to systems with a spin-flip length much larger than the lateral dimensions of the system, i.e., there are no spatial variations of accumulated spins inside the dot, though no assumptions are made about the relation between the spin-flip time $`\tau _{\mathrm{sf}}`$ and $`\tau _d`$. We stress that all the above constraints form a situation that is experimentally accessible with present techniques. Vertical structures like metallic thin-film heterostructures, epitaxial layers of half-metallic chromium dioxide or diluted magnetic semiconductors are equally suitable for the observation of the effects discussed here. In fact, a uniform spin accumulation was very recently achieved in a normal metallic island coupled to magnetic leads.zaf03

General theoretical formulations of semiclassical spin transport, based on quasiclassical Green’s functions, have been developed by Brataas, Nazarov and Bauer for normal-ferromagnetic systemsbra00  and by Belzig et al. for normal-ferromagnetic-superconducting systems.bel00  Here, we follow the approach of Samuelsson and Büttiker sam02  extended to spin-dependent transport. This simple theory is equivalent to Refs. bel00, ; bra00,  for the geometry in Fig 1(a). It is directly formulated in terms of distribution functions and allows us to treat the current and the current correlations within a unified framework. Because the proximity effect is suppressed, the electrons and holes are independent sam02  and the system is mapped onto the equivalent circuit of Fig. 1(b) in which the electron and hole subsystems (with opposite spins) are connected via Andreev reflection. Spin relaxation in the normal conductor gives rise to a spin-flip current between the spin subsystems with the same quasiparticle charge. As a result, we can describe the conductor with four nodes \[two nodes ($``$ and $``$) for each quasiparticle type (electrons and holes)\] connected to each other and to the ferromagnetic reservoirs via spin dependent resistances. We model the spin dependent conductances of the FN contacts as $`g_\alpha ^,=g(1\pm \eta _\alpha )/2`$, where $`\eta _\alpha `$ is the spin polarization of the lead $`\alpha `$ ($`\alpha =\{1,2\}`$): $`\eta _\alpha =(\nu _\alpha ^{}\nu _\alpha ^{})/(\nu _\alpha ^{}+\nu _\alpha ^{})`$, where $`\nu _\alpha ^\sigma `$ is the density of states of the ferromagnetic lead close to the Fermi energy ($`\sigma =\{,\}`$).

Conservation of time-averaged (spin-dependent) current densities must be fulfilled at each node for every energy $`E`$. Currents flowing into the normal conductor are defined positive such that the spectral currents obey

$$\underset{\alpha }{}i_{\alpha q}^{}i_S^{}+i_{,e}=0,\underset{\alpha }{}i_{\alpha q}^{}i_S^{}i_{,e}=0,$$
(1)
where $`q=(e,h)`$ labels the carrier (electron or hole) and $``$ ($`+`$) is for $`q=e`$ ($`q=h`$). The currents through the FN and NS contacts are:

$$i_{\alpha q}^\sigma =g_\alpha ^\sigma (f_{\alpha q}\widehat{f}_q^\sigma ),i_S^\sigma =g_S(\widehat{f}_e^{\overline{\sigma }}\widehat{f}_h^\sigma )/2,$$
(2)
where $`f_{\alpha q}(E)=\{1+\mathrm{exp}([EeV_\alpha ]/kT)\}^1`$ is the equilibrium distribution function at reservoir $`\alpha `$. Note that $`f_{\alpha q}`$ is spin independent as spin relaxation mechanisms are usually highly efficient in ferromagnets. In Eq. (2), $`\widehat{f}_q^\sigma (E)`$ is the nonequilibrium distribution function of the normal conductor to be calculated from Eq. (1). The spin-flip current $`i_{,q}=\gamma (\widehat{f}_q^{}\widehat{f}_q^{})`$ equilibrates the spin distributions. Here, $`\gamma =h\nu _0/\tau _{\mathrm{sf}}`$ ($`\nu _0`$ is the density of states of the conductor at the Fermi energy and $`\tau _{\mathrm{sf}}`$ the spin-flip time).bra00

## III Conductance matrix

The total current flowing out of reservoir $`\alpha `$ is

$$_\alpha =_{\alpha e}^{}+_{\alpha e}^{}_{\alpha h}^{}_{\alpha h}^{},$$
(3)
where $`_{\alpha q}^\sigma =(e/h)𝑑Ei_{\alpha q}^\sigma (E)`$, from which we determine the conductance matrix. The conductance $`G_{\alpha \beta }`$ relates the current $`_\alpha `$ with a voltage change at lead $`\beta `$. In the case $`V_1>0`$ and $`V_2=V_S=0`$, the interesting conductance is $`G_{21}=(h/2e^2)(_2/V_1)`$:

$$G_{21}=\frac{g^2[g(\eta _{1}^{}{}_{}{}^{2}+\eta _{2}^{}{}_{}{}^{2}2)4\gamma 2\eta _1\eta _2g_S]}{g^2[4(\eta _1+\eta _2)^2]+8\gamma g_S+4g\left(2\gamma +g_S\right)}.$$
(4)
In a normal system, $`G_{21}`$ is always negative.but86  On the contrary, in our system we can have $`G_{21}>0`$, which is the manifestation of the Andreev drag effect.pha02  Figure 2 shows the behavior of $`G_{21}`$ as a function of $`2g/g_S`$. In the parallel (P) case, $`\eta _1=\eta _2=\eta `$ and $`G_{21}`$ is negative \[Fig. 2(a)\]. The maximum (absolute) value of $`G_{21}`$ is attained at $`\eta =1`$ and decreases with increasing values of $`\gamma `$ \[Fig. 2(b)\]. In the antiparallel (AP) configuration ($`\eta _1=\eta _2=\eta `$), a *reversal* of the current at F2 occurs for low values of $`g/g_S`$ and $`\gamma `$ \[see Fig. 2(c) and Fig. 2(d)\]. Note that the drag effect can result already for a small value of $`\eta `$. In fact, crossed Andreev reflections start to dominate when $`\eta ^2>(g/g_S)/(1+g/g_S)`$ at $`\gamma =0`$. In particular, for the fully polarized case one always finds $`G_{21}>0`$ *regardless of the ratio $`g/g_S`$*. When $`g_Sg`$, we find that $`G_{21}>0`$ holds for any value of $`\eta >0`$, in agreement with Ref. jed99, . Increasing $`\gamma `$, the spin relaxation processes randomize the spin distributions and eventually lead to $`G_{21}<0`$, \[the curves tend to the $`\eta =0`$ case of Fig. 2(b)\]. The condition for the drag effect is:

$$2\gamma <\eta ^2g_Sg(1\eta ^2).$$
(5)
This is a central result of our work.

Interestingly, the Andreev drag effect implies that even the equilibrium Nyquist-Johnson noise yields *positive* current cross-correlations $`P_{21}=4k_BTG_{21}`$ in contrast to normal conductors. Nevertheless, energy dissipation is ensured as the eigenvalues of the conductance matrix are semipositive.

Another striking result is that for $`\eta =1`$ there are two linearly independent voltage configurations that yield *zero* eigenvalue of the conductance matrix. One is the gauge solution ($`V_1=V_2=V_S`$). The nontrivial arrangement is $`V_1=V_2=V`$ and $`V_S=0`$. Here, F1 injects electrons whereas an equivalent flow of incident holes stems from F2. Complete Andreev reflection yields then zero net current.

## IV Magnetization

Out of equilibrium for $`V_1>0`$ and $`V_2=V_S=0`$, an accumulation of majority spins forms within the conductor. Before applying the voltage, the conductor has zero magnetization. In the AP case, the polarization excess is given by

$$M=\frac{𝑑E_q(\widehat{f}_q^{}\widehat{f}_q^{})}{𝑑E_q(\widehat{f}_q^{}+\widehat{f}_q^{})}=\frac{g\eta }{2\gamma +g+g_S}.$$
(6)
It approaches $`\eta `$ for $`gg_S,\gamma `$ and vanishes for $`\gamma g`$ or $`g_Sg`$.

## V Cross correlations

To obtain additional insight into the drag effect, we investigate the shot noise of the current. Of particular interest are the correlations between currents flowing in the two ferromagnetic reservoirs. In purely normal systems, the current cross-correlations are manifestly negative.but92  However, in hybrid superconducting structures it was shown that Andreev reflections can cause a reversal of the sign of the cross-correlations. Poscorr  The positive correlations in few-mode conductors were also shown to be enhanced by ferromagnetic leads.tad02  Here, we focus on the cross correlations in the AP configuration, which is the most interesting case.

The contacts act as generators of fluctuations $`\delta i_{\alpha q}^\sigma `$ and $`\delta i_S^\sigma `$. As a result, the spin-dependent distribution function of the normal conductor fluctuates itself. The total fluctuating current is

$$\mathrm{\Delta }i_{\alpha q}^\sigma =\delta i_{\alpha q}^\sigma g_\alpha ^\sigma \delta \widehat{f}_q^\sigma ,\mathrm{\Delta }i_S^\sigma =\delta i_S^\sigma +g_S(\delta \widehat{f}_e^{\overline{\sigma }}\delta \widehat{f}_h^\sigma ).$$
(7)
We also take into account the fluctuations $`\delta i_{,q}`$ due to spin relaxation, giving the fluctuating current

$$\mathrm{\Delta }i_{,q}=\delta i_{,q}+\gamma (\delta \widehat{f}_q^{}\delta \widehat{f}_q^{}).$$
(8)
The next step is to determine the fluctuations of the distribution functions of the conductor, $`\delta \widehat{f}_q^\sigma `$. In the zero frequency limit, the fluctuating currents are conserved ($``$ for $`q=e,h`$):

$`{\displaystyle \underset{\alpha }{}}\mathrm{\Delta }i_{\alpha q}^{}\mathrm{\Delta }i_S^{}+\mathrm{\Delta }i_{,q}=0,`$ (9a)
$`{\displaystyle \underset{\alpha }{}}\mathrm{\Delta }i_{\alpha q}^{}\mathrm{\Delta }i_S^{}\mathrm{\Delta }i_{,q}=0.`$ (9b)
The fluctuations of the charge current at contact $`\alpha `$ are found from $`\mathrm{\Delta }i_\alpha =\mathrm{\Delta }i_{\alpha e}^{}+\mathrm{\Delta }i_{\alpha e}^{}\mathrm{\Delta }i_{\alpha h}^{}\mathrm{\Delta }i_{\alpha h}^{}`$.

The noise power of the current cross-correlations is

$$P_{21}=2𝑑t\mathrm{\Delta }_2(t)\mathrm{\Delta }_1(0).$$
(10)
The contacts emit fluctuations independently so that

$$\delta i_{\alpha q}^\sigma (E,t)\delta i_{\alpha ^{}q^{}}^\sigma ^{}(E^{},t^{})=\frac{h}{e}\overline{\delta }S_{\alpha q}^\sigma ,$$
(11)
with $`\overline{\delta }\delta _{\alpha \alpha ^{}}\delta _{qq^{}}\delta _{\sigma \sigma ^{}}\delta (EE^{})\delta (tt^{})`$. Here, $`S_{\alpha q}^\sigma `$ is the fluctuation power:sam02 ; but92

$`S_{\alpha q}^\sigma =eg_\alpha ^\sigma [f_{\alpha q}^\sigma +\widehat{f}_q^\sigma 2\widehat{f}_q^\sigma f_{\alpha q}^\sigma \mathrm{\Gamma }(f_{\alpha q}^\sigma \widehat{f}_q^\sigma )^2],`$ (12a)
$`S_S^\sigma =eg_S[\widehat{f}_e^\sigma +\widehat{f}_h^{\overline{\sigma }}2\widehat{f}_h^{\overline{\sigma }}\widehat{f}_e^\sigma R(\widehat{f}_e^\sigma \widehat{f}_h^{\overline{\sigma }})^2]/2,`$ (12b)
$`S_q=e\gamma [\widehat{f}_q^{}+\widehat{f}_q^{}2\widehat{f}_q^{}\widehat{f}_q^{}],`$ (12c)
where we have assumed that the statistics of spin-flip events is Poissonian.mis03

The above system of equations is complete and can be solved for any set of applied voltages, polarization alignments and spin-flip rates. We calculate $`P_{21}`$ for the bias configuration given by $`V_1>0`$, $`V_2=V_S=0`$ in the AP case at $`k_BT=0`$. We present only the analytical result in limiting cases and show the general result in Fig. 3. In Fig. 3(a), we plot $`P_{21}`$ as a function of the conductance ratio $`2g/g_S`$ for various polarizations $`\eta `$, $`\gamma =0`$ and $`\mathrm{\Gamma }=R=1`$. (Qualitatively, the results are independent of $`\mathrm{\Gamma }`$ and $`R`$.) For $`\eta =0`$, the cross-correlations are negative. As $`\eta `$ is turned on, $`P_{21}`$ can become positive. With increasing $`g/g_S`$, $`P_{21}`$ exhibits a crossover from negative to positive values at $`2g/g_S1`$. Indeed, for dominating NS conductance (small $`g/g_S`$) we find

$$P_{21}=p_0\eta ^2[\mathrm{\Gamma }+2\eta ^2(3\mathrm{\Gamma }2)]/16+𝒪(g/g_S),$$
(13)
which is manifestly positive for $`\eta >0`$ ($`p_04gV_1e^3/h`$). In this case, the electron and hole distributions within the conductor are almost identical.

On the other hand, for dominating FN conductance (small $`g_S/g`$) the cross-correlations become

$$P_{21}=p_0(1\eta ^2)[\mathrm{\Gamma }2+\eta ^2(3\mathrm{\Gamma }2)]/8+𝒪(g_S/g),$$
(14)
which are negative for $`\eta <1`$. Notice that for $`\eta =0`$ this expression yields the shot noise for a chaotic dot coupled to two normal leads. In the fully polarized case, the cross-correlations are *always* positive independently of $`g/g_S`$.

For illustrative purposes we show in the inset of Fig. 3(a) the “phase diagram” $`(\eta ,g/g_S)`$. The solid line marks the crossover between positive and negative cross-correlations. With increasing $`\eta `$, the crossover point shifts toward larger values of $`2g/g_S`$. For comparison, we also plot a dashed curve that marks the transition to the Andreev drag effect. Notice that the conditions $`P_{21}>0`$ and $`G_{21}>0`$ are closely related, but not identical. This is the case since the scattering processes in general contribute in different ways for the current and noise.sam02

A finite source of spin relaxation leads to the suppression of positive $`P_{21}`$. To see this, we show in Fig. 3(b) the cross-correlations for the fully polarized case for different values of $`\gamma `$. For nonzero values of $`\gamma `$, $`P_{21}`$ becomes negative with increasing $`2g/g_S`$. Eventually, when $`\gamma \mathrm{}`$, we recover the unpolarized case ($`\eta =0`$) of Fig. 3(a).

## VI Conclusions

We have investigated the transport properties of a normal conductor attached to ferromagnetic leads and coupled to a superconductor in the subgap regime. The proposed geometry, which includes generic elastic and spin-flip scattering, is shown to exhibit a pronounced Andreev drag effect. We have demonstrated that the equilibrium and transport current correlations are strongly sensitive to crossed Andreev reflections. The theory presented here applies to a wide range of geometries and should therefore advance experimental observation of the Andreev drag effect.

Note added in proof. Recently, we became aware of the work of Lambert et al.lam03  who compute numerically the conductances for a metallic diffusive geometry similar to Fig. 1. They also conclude that a normal conductor inserted between the ferromagnetic and the superconducting contacts leads to an enhancement of the Andreev drag effect. In contrast to our work, the role of spin-flip, the effect of differing lead polarizations, and the shot noise are not investigated.

## Acknowledgments

We thank V. Chandrasekhar and B.J. van Wees for helpful discussions. We acknowledge support from the Swiss NSF through the program MANEP. D.S. and R.L. were also supported by the EU RTN under Contract No. HPRN-CT-2000-00144, Nanoscale Dynamics.