text:
# 300x Faster Matlab using MatlabMPI This work is sponsored by the High Performance Computing Modernization Office, under Air Force Contract F19628-00-C-0002. Opinions, interpretations, conclusions and recommendations are those of the author and are not necessarily endorsed by the Department of Defense.

## 1 Introduction

Matlab  is the dominant programming language for implementing numerical computations and is widely used for algorithm development, simulation, data reduction, testing and system evaluation. The popularity of Matlab is driven by the high productivity that is achieved by users because one line of Matlab code can typically replace ten lines of C or Fortran code. Many Matlab programs can benefit from faster execution on a parallel computer, but achieving this goal has been a significant challenge. There have been many previous attempts to provide an efficient mechanism for running Matlab programs on parallel computers . These efforts of have faced numerous challenges and none have received widespread acceptance.

In the world of parallel computing the Message Passing Interface (MPI)  is the de facto standard for implementing programs on multiple processors. MPI defines C and Fortran language functions for doing point-to-point communication in a parallel program. MPI has proven to be an effective model for implementing parallel programs and is used by many of the worlds’ most demanding applications (weather modeling, weapons simulation, aircraft design, and signal processing simulation).

MatlabMPI consists of a set of Matlab scripts that implements a subset of MPI and allows any Matlab program to be run on a parallel computer. The key innovation of MatlabMPI is that it implements the widely used MPI “look and feel” on top of standard Matlab file I/O, resulting in a “pure” Matlab implementation that is exceedingly small ($``$250 lines of code). Thus, MatlabMPI will run on any combination of computers that Matlab supports.

The next section describes the implementation and functionality provided by MatlabMPI. Section three presents results on the bandwidth performance of the library from several parallel computers. Section four uses an image processing application to show the scaling performance and the very high software price performance that can be achieved using MatlabMPI. Section five presents our conclusions and plans for future work.

## 2 Implementation

The central innovation of MatlabMPI is the use of Matlab’s built in file I/O capabilities, which allow any Matlab variable (matrices, arrays, structures, …) to be written and read by Matlab running on any machine, and eliminates the need to write complicated buffer packing and unpacking routines. The approach used in MatlabMPI is extremely simple and is illustrated in Figure 1. The sender saves a Matlab variable to a data file and when the file is complete the sender touches a lock file. The receiver continuously checks for the existence of the lock file; when it is exists the receiver reads in the data file and then deletes both the data file and the lock file. An example of a basic send and receive Matlab program is shown below

> ```
> MPI_Init;                           % Initialize MPI.
> comm = MPI_COMM_WORLD;              % Create communicator.
> comm_size = MPI_Comm_size(comm);    % Get size.
> my_rank = MPI_Comm_rank(comm);      % Get rank.
> source = 0;                         % Set source.
> dest = 1;                           % Set destination.
> tag = 1;                            % Set message tag.
> if(comm_size == 2)                  % Check size.
>   if (my_rank == source)            % If source.
>     data = 1:10;                    % Create data.
>     MPI_Send(data,dest,tag,comm);   % Send data.
>   end
>   if (my_rank == dest)              % If destination.
>     data=MPI_Recv(source,tag,comm); % Receive data.
>   end
> end
> MPI_Finalize;                       % Finalize Matlab MPI.
> exit;                               % Exit Matlab
>
> ```

The structure of the above program is very similar to MPI programs written in C or fortran, but with the convenience of a high level language. The first part of the program sets up the MPI world; the second part differentiates the program according to rank and executes the communication; the third part closes down the MPI world and exits Matlab. If the above program were written in a matlab file called SendReceive.m, it would be executed by calling the following command from the Matlab prompt:

```
    eval( MPI_Run(’SendReceive’,2,machines) );

```

Where the machines argument can be any of the following forms:

> Run on the local host.
>
> Run on a multi processors.
>
> Run on a multi processors and communicate using dir1 and dir2, which must be visible to both machines.

The MPI\_Run command launches a Matlab script on the specified machines with output redirected to a file. If the rank=0 process is being run on the local machine, MPI\_Run returns a string containing the commands to initialize MatlabMPI, which allows MatlabMPI to be invoked deep inside of a Matlab program in a manner similar to fork-join model employed in OpenMP.

The SendReceive example illustrates the basic six MPI functions (plus MPI\_Run) that have been implemented in MatlabMPI

> Runs a matlab script in parallel.
>
> Initializes MPI.
>
> Gets the number of processors in a communicator.
>
> Gets the rank of current processor within a communicator.
>
> Sends a message to a processor (non-blocking).
>
> Receives message from a processor (blocking).
>
> Finalizes MPI.

For convenience, three additional MPI functions have also been implemented along with two utility functions that provide important MatlabMPI functionality, but are outside the MPI specification

> Function to kill all matlab jobs started by MatlabMPI.
>
> Broadcast a message (blocking).
>
> Returns a list of all incoming messages.
>
> MatlabMPI function to prevent messages from being deleted (useful for debugging).
>
> MatlabMPI function to delete all files created by MatlabMPI.

MatlabMPI handles errors the same as Matlab, however running hundreds of copies does bring up some additional issues. If an error is encountered and the Matlab script has an “exit” statement then all the Matlab processes will die gracefully. If a Matlab job is waiting for a message that will never arrives than it needs to be killed with the MPI\_Abort command. In this situation, MatlabMPI can leave files which need to be deleted with the MatMPI\_Delete\_all command.

On shared memory systems, MatlabMPI only requires a single Matlab license since any user is allowed to launch many Matlab sessions. On a distributed memory system, MatlabMPI requires one Matlab license per machine. Because MatlabMPI uses file I/O for communication, there must also be a directory that is visible to every machine (this is usually also required in order to install Matlab). This directory defaults to the directory that the program is launched from.

## 3 Bandwidth

The vast majority of potential Matlab applications are “embarrassingly” parallel and require minimal performance out of MatlabMPI. These applications exploit coarse grain parallelism and communicate rarely (if at all). Never-the-less, measuring the communication performance is useful for determining which applications are most suitable for MatlabMPI.

MatlabMPI has been run on several Unix platforms. It has been benchmarked and compared to the performance of C MPI on the SGI Origin2000 housed at Boston University. These results indicate that for large messages ($``$1 MByte) MatlabMPI is able to match the performance of C MPI (see Figure 2). For smaller messages, MatlabMPI is dominated by its latency ($``$35 milliseconds), which is significantly larger than C MPI. These results have been reproduced using a SGI Origin2000 and a Hewlett Packard workstation cluster (connected with 100 Mbit ethernet) at Ohio State University (see Figure 3).

The above bandwidth results were all obtained using two processors engaging in bi-directional sends and receives. Such a test does a good job of testing the individual links on an a multi-processor system. To more broadly test the interconnect the send receive benchmark is run on an eight node (16 cpu) Linux cluster connected with Gigabit ethernet (see Figure 4). These results are shown for one and 16 cpus. MatlabMPI is able to maintain high bandwidth even when multiple processors are communicating by allowing each processor to have its own receive directory. By cross mounting all the disks in a cluster each node only sees the traffic directed to it, which allows communication contention to be kept to a minimum.

## 4 Scaling

To further test MatlabMPI a simple image filtering application was written. This application abstracts the key computations that are used in many of our DoD sensor processing applications (e.g. wide area Synthetic Aperture Radar). The application executed repeated 2D convolution on a large image (1024 x 128,000 $``$2 GBytes). This application was run on a large shared memory parallel computer (the SGI Origin2000 at Boston University) and achieved speedups greater than 64 on 64 processors; showing the classic super-linear speedup (due to better cache usage) that comes from breaking a very large memory problem into many smaller problems (see Figure 5).

To further test the scalability, the image processing application was run with a constant load per processor (1024 x 1024 image per processor) on a large shared/distributed memory system (the IBM SP2 at the Maui High Performance Computing Center). In this test, the application achieved a speedup of $``$300 on 304 CPUs as well achieving $``$15% of the theoretical peak (450 Gigaflops) of the system (see Figure 6).

The ultimate goal of running Matlab on parallel computers is to increase programmer productivity and decrease the large software cost of using HPC systems. Figure 7) plots the software cost (measured in Software Lines of Code or SLOCs) as a function of the maximum achieved performance (measured in units of single processor peak) for the same image filtering application implemented using several different libraries and languages (VSIPL, MPI, OpenMP, using C++, C, and Matlab ). These data show that higher level languages require fewer lines to implement the same level of functionality. Obtaining increased peak performance (i.e. exploiting more parallelism) requires more lines of code. MatlabMPI is unique in that it achieves a high peak performance using a small number of lines of code.

Two useful metrics we have developed for measuring software productivity are Gigaflops/SLOC and CPUs/SLOC. Gigaflops/SLOC will increase in a given application if it can be made to run faster or can be implemented with fewer lines of code. CPUs/SLOC will increase if the number of processors that can be effectively used in an application can be increased without increasing the number of lines of code. The test application does extremely well in both of these measures, achieving 0.85 Gigaflops/SLOC and 4.4 CPUs/SLOC, both of which are the highest values recorded for any application.

## 5 Conclusions and Future Work

The use of file I/O as a parallel communication mechanism is not new and is now increasingly feasible with the availability of low cost high speed disks. The extreme example of this approach are the now popular Storage Area Networks (SAN), which combine high speed routers and disks to provide server solutions. Although using file I/O increases the latency of messages it normally will not effect the bandwidth. Furthermore, the use of file I/O has several additional functional advantages which make it easy to implement large buffer sizes, recordable messages, multi-casting, and one-sided messaging. Finally, the MatlabMPI approach is readily applied to any language (e.g. IDL, Python, and Perl).

MatlabMPI demonstrates that the standard approach to writing parallel programs in C and Fortran (i.e. using MPI) is also valid in Matlab. In addition, by using Matlab file I/O, it was possible to implement MatlabMPI entirely within the Matlab environment, making it instantly portable to all computers that Matlab runs on. Most potential parallel Matlab applications are trivially parallel and don’t require high performance. Never-the-less, MatlabMPI can match C MPI performance on large messages. The simplicity and performance of MatlabMPI makes it a very reasonable choice for programmers that want to speed up their Matlab code on a parallel computer.

MatlabMPI provides the highest productivity parallel computing environment available. However, because it is a point-to-point messaging library, a significant amount code of must be added to any application in order to do basic parallel operations. In the test application presented here, the number of lines of Matlab code increased from 35 to 70. While a 70 line parallel program is extremely small, it represents a significant increase over the single processor case. Our future work will aim at creating higher level objects (e.g. distributed matrices) that will eliminate this parallel coding overhead (see Figure 8). The resulting “Parallel Matlab Toolbox” will be built on top of the MatlabMPI communication layer, and will allow a user to achieve good parallel performance without increasing the number lines of code.

## Acknowledgments

The authors would like to thank the sponsorship of the DoD High Performance Computing Modernization Office, and the staff at the supercomputing centers at Boston University, Ohio State University and the Maui High Performance Computing Center.

## Appendix A Parallel Image Filtering Test Application

```
MPI_Init;  % Initialize MPI.
comm = MPI_COMM_WORLD;  % Create communicator.
comm_size = MPI_Comm_size(comm);  % Get size.
my_rank = MPI_Comm_rank(comm);  % Get rank.
% Do a synchronized start.
starter_rank = 0;
delay = 30;  % Seconds
synch_start(comm,starter_rank,delay);
n_image_x = 2.^(10+1)*comm_size;  % Set image size (use powers of 2).
n_image_y = 2.^10;
n_point = 100;  % Number of points to put in each sub-image.
% Set filter size (use powers of 2).
n_filter_x = 2.^5;
n_filter_y = 2.^5;
n_trial = 2;  % Set the number of times to filter.
% Computer number of operations.
total_ops = 2.*n_trial*n_filter_x*n_filter_y*n_image_x*n_image_y;
if(rem(n_image_x,comm_size) ~= 0)
 disp(’ERROR: processors need to evenly divide image’);
 exit;
end
disp([’my_rank: ’,num2str(my_rank)]);  % Print rank.
left = my_rank - 1;  % Set who is source and who is destination.
if (left < 0)
  left = comm_size - 1;
end
right = my_rank + 1;
if (right >= comm_size)
  right = 0;
end
tag = 1;  % Create a unique tag id for this message.
start_time = zeros(n_trial);  % Create timing matrices.
end_time = start_time;
zero_clock = clock;  % Get a zero clock.
n_sub_image_x = n_image_x./comm_size;  % Compute sub_images for each processor.
n_sub_image_y = n_image_y;
% Create starting image and working images..
sub_image0 = rand(n_sub_image_x,n_sub_image_y).^10;
sub_image = sub_image0;
work_image = zeros(n_sub_image_x+n_filter_x,n_sub_image_y+n_filter_y);
% Create kernel.
x_shape = sin(pi.*(0:(n_filter_x-1))./(n_filter_x-1)).^2;
y_shape = sin(pi.*(0:(n_filter_y-1))./(n_filter_y-1)).^2;
kernel = x_shape.’ * y_shape;
% Create box indices.
lboxw = [1,n_filter_x/2,1,n_sub_image_y];
cboxw = [n_filter_x/2+1,n_filter_x/2+n_sub_image_x,1,n_sub_image_y];
rboxw = [n_filter_x/2+n_sub_image_x+1,n_sub_image_x+n_filter_x,1,n_sub_image_y];
lboxi = [1,n_filter_x/2,1,n_sub_image_y];
rboxi = [n_sub_image_x-n_filter_x/2+1,n_sub_image_x,1,n_sub_image_y];
start_time = etime(clock,zero_clock);  % Set start time.
% Loop over each trial.
for i_trial = 1:n_trial
  % Copy center sub_image into work_image.
  work_image(cboxw(1):cboxw(2),cboxw(3):cboxw(4)) = sub_image;
  if (comm_size > 1)
    ltag = 2.*i_trial;    % Create message tag.
    rtag = 2.*i_trial+1;
    % Send left sub-image.
    l_sub_image = sub_image(lboxi(1):lboxi(2),lboxi(3):lboxi(4));
    MPI_Send(  left, ltag, comm, l_sub_image );
    % Receive right padding.
    r_pad = MPI_Recv( right, ltag, comm );
    work_image(rboxw(1):rboxw(2),rboxw(3):rboxw(4)) = r_pad;
    % Send right sub-image.
    r_sub_image = sub_image(rboxi(1):rboxi(2),rboxi(3):rboxi(4));
    MPI_Send( right, rtag, comm, r_sub_image );
    % Receive left padding.
    l_pad = MPI_Recv( left, rtag, comm );
    work_image(lboxw(1):lboxw(2),lboxw(3):lboxw(4)) = l_pad;
  end
  work_image = conv2(work_image,kernel,’same’);   % Compute convolution.
  % Extract sub_image.
  sub_image = work_image(cboxw(1):cboxw(2),cboxw(3):cboxw(4));
end
end_time = etime(clock,zero_clock);  % Get end time for the this message.
total_time = end_time - start_time  % Print the results.
% Print compute performance.
gigaflops = total_ops / total_time / 1.e9;
disp([’GigaFlops: ’,num2str(gigaflops)]);
MPI_Finalize;  % Finalize Matlab MPI.
exit;

```


format: text

text:
# Spontaneous Fluxon Formation in Annular Josephson Tunnel JunctionsSubmitted to Phys. Rev. B.

## I Introduction

As the early universe cooled it is believed to have undergone a series of spontaneous phase transitions, whose inhomogeneities (monopoles, cosmic strings, domain walls) have observable consequences, for structure formation in particular. These defects appear because the correlation length $`\xi `$ of the field (or fields) whose expectation value is the order parameter is necessarily finite for a transition that is implemented in a finite time, whether it be continuous or not.

It is difficult to determine the microscopic dynamics of such fields but, using only simple causal arguments Kibble made estimates of this early field ordering, and the density of topological defects produced at Grand Unified Theory (GUT) transitions at $`10^{35}s`$. Unfortunately, because the nature of the field theories is not known with any reliability, and the effects of their evolution are not visible until the decoupling of the radiation and matter $`10^6yrs`$ later, it is impossible to provide unambiguous checks of these predictions. However, causality is such a fundamental notion that Zurek suggested that identical causal arguments, with similar predictions, were applicable to condensed matter systems for which direct experiments on defects could be performed. In addition to their intrinsic interest for a better understanding of the dynamics of transitions in condensed matter, the hope is that successful tests of these predictions can lead to a better understanding of phase transitions in quantum fields.

Several experiments in condensed matter systems have already been performed to test the Zurek-Kibble predictions, with mixed results. It is these predictions that we have tested here, using annular Josephson tunnel junctions, for which the defects are fluxons.

This paper is organised as follows. In the next section we discuss the Zurek-Kibble (ZK) scenario. Successive sections then give the predictions for annular Josephson Tunnel Junctions (AJTJs), show how the fluxons are measured, the nature of the samples, describe the experimental setup and, finally, present the measurements and their agreement with the predictions. An abbreviated description of the experiment has been given elsewhere. Details of the theoretical analysis can be found in , in which an earlier experiment by two of us (R.M and J.M) was analysed, to demonstrate its compatibility with the ZK analysis, even though it had not been performed with this in mind.

## II Zurek-Kibble Causality

Consider a system with critical temperature $`T_c`$, cooled through that temperature so that, if $`T(t)`$ is the temperature at time $`t`$, then $`T(0)=T_c`$. $`\dot{T}(0)=T_c/\tau _Q`$ defines the quench time $`\tau _Q`$.

There are several ways of formulating the Zurek-Kibble causality bounds, but they all depend on the fact that, as the transition begins to be implemented, there is a maximum speed $`c(t)=c(T(t))`$ at which the system can become ordered. For relativistic quantum field theory (QFT) $`c`$ is the speed of light. For superfluids, $`c(t)`$ is the speed of second sound, vanishing at $`t=0`$. For Josephson tunnel Junctions $`c(t)`$, which depends on the nature of the junction, is the Swihart velocity.

Suppose that the equilibrium (adiabatic) correlation length $`\xi _{ad}(t)=\xi _{ad}(T(t))`$ diverges near $`t=0`$ as

$`\xi _{ad}(t)=\xi _0\left|{\displaystyle \frac{t}{\tau _Q}}\right|^\nu .`$

However, the true non-equilibrium correlation length $`\xi (t)`$ can only change so much in a finite time, and does not diverge. Kibble and Zurek made two assumptions:

* Firstly, the correlation length $`\overline{\xi }`$ of the fields that characterizes the onset of order is the equilibrium correlation length $`\overline{\xi }=\xi _{ad}(\overline{t})`$ at some appropriate time $`\overline{t}`$.
* Secondly, we can measure $`\overline{\xi }`$ experimentally by measuring the number of defects, assuming that the defect separation $`\xi _{def}=O(\overline{\xi })`$.

There are several ways to estimate $`\overline{t}`$ (the ’causal time’), explicit or implicit in the early work of Zurek. Most simply,

* $`\xi (t)`$ cannot grow faster than $`c(t)`$. This is true both before and after the transition. That is, $`\overline{t}`$ is defined by the condition that $`\dot{\xi }_{ad}(\overline{t})c(\overline{t})`$.
* The relaxation time for the long wavelength modes $`\tau (t)`$ is defined by $`c(t)=\xi _{ad}(t)/\tau (t)`$. From this viewpoint $`\overline{t}`$ is that time when we can return to an adiabatic regime, $`\overline{t}\tau (\overline{t})`$.

In simple systems these estimates agree, up to numerical factors approximately unity. They give $`\overline{t}`$ of the form $`\overline{t}\tau _Q^{1\gamma }\tau _0^\gamma `$, where $`\tau _0\tau _Q`$ is the cold relaxation time of the longest wavelength modes, and the critical exponent $`\gamma `$ depends upon the system. As a result, $`\tau _Q\overline{t}\tau _0`$. Identifying the initial domain size and defect separation as in the second ZK assumption then gives

$$\overline{\xi }\xi _{ad}(\overline{t})=\xi _0\left(\frac{\tau _Q}{\tau _0}\right)^\sigma \xi _0,$$
(1)
where $`\sigma =\gamma \nu `$. This is very large on the scale of cold defects which shrink to size $`\xi _{ad}(T_{fin})=O(\xi _0)`$, where $`T_{fin}`$ is the final temperature, and $`\xi _0`$ is determined from the microscopic dynamics. We term $`\sigma `$ the Zurek-Kibble (ZK) scaling index.

The analysis above is for large systems, of linear size $`L\overline{\xi }`$. For the experiment that we shall describe below, of Annular Josephson Tunnel Junctions (AJTJs) of circumference $`C`$, we find that is $`C<\overline{\xi }`$. In fact, we expect the conclusions to be equally valid for small systems, for which the relevant time might, incorrectly, seem to be the time when the coherence length becomes smaller than the system. The reason is that the causal bounds (1) are to be thought of as a short-hand for the underlying dynamics. At the microscopic level, causality along the lines above is not explicit, although encoded in the relevant dynamical equations. The picture is rather one of order being established through the growth of the amplitudes of long-wavelength instabilities. The earliest time at which we can identify defects from this viewpoint is when the order parameters have achieved their equilibrium magnitudes. Qualitatively, for simple models this time is in good agreement with the causal time $`\overline{t}`$ above. There is no real surprise in this. It has been shown by one of us that, in general, the causal time and distance scales $`\overline{t}`$ and $`\overline{\xi }=\xi (\overline{t})`$ are just as we would expect from dimensional analysis (in the mean-field approximation) and unstable modes grow exponentially, whereby the dependence of the causal time (and corresponding defect density) on the microscopic parameters is only logarithmic. In the same way, the distance between field zeroes has the same scaling dependence on $`\tau _Q`$ as $`\overline{\xi }`$ of (1), up to logarithms.

As for the production of defects, they are transitional regions between different system ground-states. For superfluids like $`{}_{}{}^{4}He`$ and superconductors these transitional regions are flagged by zeroes of the scalar order-parameter field. For the case of JTJs we shall see that this is generalised to zeroes (mod $`2\pi `$) of the order-parameter field. Any field crossing zero (mod $`2\pi `$) has the potential to mature into a defect. However, only when the transition is complete will the field configurations in the vicinity of the zeroes have the energy profile of a classical defect, the solution to the classical field equations. Thus, before the causal time we now have a picture in which there is a fractal thermal fuzz of potential defects, whose density depends on the scale at which we look. By the causal time some of these have developed into the (scale-independent) defects that we see subsequently (see ). Because, from this viewpoint, we have many proto-defects jockeying to become the real thing, the relevant scale to compare to the system size is not $`\overline{\xi }`$ but $`\xi _0`$, as before. Equivalently, counting field zeroes (mod $`2\pi `$, or not) depends on the short-distance behaviour of the field correlations. We do not expect a problem as long as $`C\xi _0`$, as is always the case.

We conclude by noting that this is a very different picture from that of domains freezing in as the transition is approached from above ($`t<0`$), which is how causal bounds were originally posed. What matters is that all these causal descriptions give results with the correct engineering dimensions. With this in mind, we keep the causal bound (1) as a convenient mnemonic.

## III The ZK Predictions For Fluxons

The order parameter of a Josephson tunnel junction at temperature $`T<T_c`$ is the phase difference $`\varphi `$ of the macroscopic superconducting quantum mechanical wave functions across the barrier. Using a Lagrangian formalism Grønbech-Jensen et al.  have shown that, for an annular JTJ with a distributed bias current $`\gamma `$, $`\varphi `$ obeys to the following perturbed Sine-Gordon equation (PSGE):

$`{\displaystyle \frac{^2\varphi }{x^2}}`$ $``$ $`{\displaystyle \frac{1}{c^2(T)}}{\displaystyle \frac{^2\varphi }{t^2}}{\displaystyle \frac{1}{\lambda _J^2(T)}}\mathrm{sin}\varphi `$ (2)
$`=`$ $`\gamma +{\displaystyle \frac{\alpha }{c^2(T)}}{\displaystyle \frac{\varphi }{t}}\beta {\displaystyle \frac{^3\varphi }{x^2t}}`$ (3)
provided the width $`\mathrm{\Delta }r`$ of the annulus, of radius $`r`$, satisfies $`\mathrm{\Delta }rr`$ and $`\mathrm{\Delta }r\lambda _J(T)`$, the Josephson coherence length. In this case $`x`$ measures the distance along the annulus, and $`c(T)`$ is the Swihart velocity. $`\alpha `$ and $`\beta `$ are the coefficients of the losses due to the tunnelling current and to the surface impedance, respectively.

The boundary conditions for Eq.(3) are periodic  and derive from fluxoid quantization : $`\varphi (x+C)=\varphi (x)+2\pi n,`$ where $`C=2\pi r`$ is the circumference of the junction and the winding number $`n`$ is an integer corresponding to the algebraic sum of fluxons trapped in the junction barrier at the normal-superconducting (N-S) transition; $`n`$ is a topological system constant, that is, only fluxon-antifluxon ($`F\overline{F}`$) pairs can be created or annihilated as long as the junction remains in the superconducting state.

The classical fluxons are the ’kinks’ of the Sine-Gordon theory. As with other models of defect formation, equation (3) is only valid once the transition is complete: therefore, we shall not use it to study the appearance of fluxons. However, it is sufficient to enable us, in the spirit of the Zurek-Kibble scenario, to identify $`\lambda _J(T)`$, diverging at $`T_c`$, as the equilibrium correlation length $`\xi _{ad}(T)`$ to be constrained by causality. Further, the Swihart velocity $`c(T)`$ (with critical slowing down at $`T=T_c`$), measures the maximum speed at which the order parameter can change.

A detailed discussion of the ZK bounds has been given elsewhere by us, and we refer the reader to these articles for more details. The JTJs in our experiment are symmetric, by which is meant that the electrodes are made of identical superconducting material with the same energy gaps and the same $`T_c`$. For such JTJs  $`\gamma =\nu =1/2`$. Therefore, at the time of their formation

$`\overline{t}=\sqrt{\tau _Q\tau _0},`$

the separation of fluxons is expected, in the ZK picture, to be

$$\overline{\xi }\xi _0\left(\frac{\tau _Q}{\tau _0}\right)^{1/4},$$
(4)
where $`\xi _0`$ and $`\tau _0`$ have the same meaning as before. It is the prediction (4) that will be tested in our experiment.

In terms of the parameters of the JTJs, $`\xi _0`$ has been inferred as

$`\xi _0=\sqrt{{\displaystyle \frac{\mathrm{}}{2e\mu _0d_saJ_c(0)}}}.`$

where $`J_c(T)`$ is the Josephson current density at temperature $`T`$. The parameter $`a`$ is given in terms of the superconductor gap energy and critical temperature and has a value between $`3`$ and $`5`$. If the thickness of the two superconducting electrodes differs, the effective thickness $`d_s`$ is the harmonic mean of the individual thicknesses. As for $`\tau _0`$, it is given as $`\tau _0=\xi _0/c_0`$, where $`c_0`$ defines the behavior $`c(t)=c_0(t/\tau _Q)^{1/2}`$of the Swihart velocity for the system near $`T=T_c`$.

## IV Measuring Fluxons

Once fluxons have appeared Eq.3 is relevant. A consequence of the periodic boundary conditions in AJTJs is that fluxons behave as relativistic particles on an infinite lossy line. In the absence of any current through the barrier and/or externally applied magnetic field the fluxons experience a flat potential and therefore are in indifferent equilibrium as far as the barrier is homogeneous and pin-hole free; in the reality, the barrier defects act as small pinning or repulsive potentials for the fluxons. Due to the losses, after the transient regime is over the fluxons are still. Unfortunately static fluxons are difficult to reveal, since according to the second Josephson law, any static phase profile $`\varphi _t(x,t)=0`$ does not alter the junction zero-voltage state that is also typical of a flat profile $`\varphi (x,t)=const.`$, corresponding to the absence of any trapped fluxon.

In contrast, whenever fluxons travel around an AJTJ they leave a clear signature on the junction current-voltage characteristic (IVC) and therefore are easily detectable. In fact, as soon as a bias current is fed to the AJTJ, the fluxons move as magnetic dipoles under the action of the resulting Lorentz force. The fluxon dynamics in long JTJs is a well known topic and has received a great deal of both theoretical and experimental attention in the last few decades. If the external bias is assumed to be uniform over the junction area, then, as a result of the balance between the externally supplied power and the internally dissipated power, the fluxons move with a constant speed: the larger is the external bias the greater is the fluxon speed, but never exceeding its relativistic limit sets by the Swihart velocity. The motion direction depends both on the current sign and fluxon polarity (i.e., whether $`n=\pm 1`$) for a given bias current. Fluxons having different polarities travel in opposite directions and are likely to annihilate when they collide at low speed.

Quantitatively, if a fluxon travels around an AJTJ having a mean circumference $`C`$ with a constant speed $`v`$, then it has a angular speed $`\omega =2\pi v/C`$ and the phase $`\varphi `$ advances of $`2\pi `$ each period $`T=2\pi /\omega =C/v`$. Therefore, according to the second Josephson equation, an average voltage $`V`$ develops across the junction equal to $`V=\frac{\mathrm{\Phi }_0}{2\pi }d\varphi /dt=\mathrm{\Phi }_0/T=\mathrm{\Phi }_0v/C`$ that is proportional to the fluxon velocity. In other words, the presence of a travelling fluxon sets the junction in a finite voltage state than can be easily measured on its IVC.

By also changing the bias current through the barrier the voltage drop changes and a new branch called Zero Field Step (ZFS) appears on the junction IVC; the ZFS represents the relationship between the applied Lorentz force (proportional to the bias current) and the fluxon speed (proportional to the voltage). When $`N`$ fluxons travel around an AJTJ, the last expression is easily generalized to give a junction voltage $`V=N\mathrm{\Phi }_0C/v`$. In the last expression $`N`$ is the total number of travelling fluxons and can be larger than the winding number $`n`$ if $`F\overline{F}`$ pairs are travelling around the annulus. Therefore, we count the number of travelling fluxons by simply measuring the voltage across the AJTJ.

These properties makes AJTJs very competitive with respect to other solid state systems proposed to test the Zurek-Kibble mechanism. Our idea is to perform a large number of N-S transitions on the same AJTJ with no external current or magnetic field; at the end of each cycle the possible spontaneously generated fluxons are static. Then we supply an external current that sets the fluxons (if any) in motion around the annulus and measure the number of travelling fluxons by a careful inspection of the junction IVC. Due to the annihilation of a fluxon-antifluxon pair, this idea works well as long as the chances to spontaneously generate two fluxons are small.

Fig.1a, Fig.1b and Fig.1c represent the IVC of the same AJTJ with no fluxon trapped, with one fluxon trapped and with two fluxons trapped, respectively. We note that with no trapped fluxons the zero voltage current is very large and only $`F\overline{F}`$ current steps appear at finite voltage. In the other two cases the supercurrent is rather small (theoretically it should be vanishingly small in ideal, pinhole free barriers) and large current branches can be observed at finite voltages corresponding to the fluxons and, possibly, $`F\overline{F}`$ pairs travelling around the junction.

## V The samples

High quality $`Nb/AlAl_{ox}/Nb`$ JTJs were fabricated on $`0.5mm`$ thick silicon substrates using the trilayer technique in which the junction is realized in the window opened in a $`SiO`$ insulator layer. Details of the fabrication process can be found in Ref.. On each $`15\times 24mm^2`$ chip four JTJs were integrated, of which three ring shaped junctions having a mean circumference $`C=500\mu m`$ and a width $`\mathrm{\Delta }r=4\mu m`$ and one $`4\times 500\mu m^2`$ overlap-type linear junction. The so called ”idle region”, i.e. the overlapping of the wiring layer onto the base electrode was about $`3\mu m`$ for all the junctions. The thicknesses of the base, top and wiring layer were $`200`$, $`80`$ and $`400nm`$, respectively.

For all samples the high quality has been inferred by a measure of the I-V characteristic at $`T=4.2K`$ . In fact, the subgap current $`I_{sg}`$ at $`2mV`$ was small compared to the current rise $`\mathrm{\Delta }I_g`$ in the quasiparticle current at the gap voltage $`V_g`$, typically $`\mathrm{\Delta }I_g>35I_{sg}`$, the gap voltage was as large as $`V_g=2.76mV`$ and the maximum critical current $`I_c`$ was larger than $`0.55\mathrm{\Delta }I_g`$ for the overlap type junction. Furthermore, the application of a strong enough external magnetic field in the barrier plane completely suppressed any Josephson structures indicating the absence of electrical shorts in the barrier. It is important to mention that i) no logarithmic singularity has been observed in the IVCs at low voltages and ii) the temperature dependence of the critical current was linear as the temperature $`T`$ approached the critical temperature $`T_C`$; both these observations assure us that the junctions are symmetric, i.e. no detectable difference can be assumed between both the energy gaps $`\mathrm{\Delta }_{1,2}`$ and the critical temperatures $`T_{C1,2}`$ of the junction electrodes $`1`$ and $`2`$ in the proximity of the barrier. The maximum Josephson current density $`J_J`$ was of the order of $`1kA/cm^{2\text{ }}`$corresponding to a specific barrier normal resistance $`\rho _N`$ of about $`200\mathrm{\Omega }\mu m^2`$.

Many samples have been measured. For clarity only two will be discussed here. The geometrical and electrical (at $`4.2K`$) parameters of the two selected annular junctions on different wafers are listed in Table I. They have the same geometry (both the base and top electrode have a hole concentric to the ring,), but differ in the critical current densities, i.e., in the normalized mean circumferences $`C/\lambda _J`$. The critical current density has been calculated from the measured quasiparticle current step, $`\mathrm{\Delta }I_g`$, at the gap voltage. The values of the barrier magnetic thickness $`\mathrm{\Lambda }=180nm`$ has been used for numeric calculations. On each same chip a linear overlap-type junction with the same width, length and idle region was used in order to measure the junction Swihart velocity $`c_0`$ with a geometry in which the effects of the self field are minimized. The value of $`c_0=1.4\times 10^7m/sec`$, due to the effect of the idle region, is 1.5 times larger than that expected for a bare junction. This value of $`c_0`$ corresponds to a value of $`0.08F/m^2`$ for the barrier effective specific capacitance.

The data in Table I show that both samples are high-quality, long ($`C\lambda _J`$) annular JTJs.

| Sample | A | B |
| --- | --- | --- |
| Mean circumference $`C(\mu m)`$ | 500 | 500 |
| Width $`\mathrm{\Delta }r(\mu m)`$ | 4 | 4 |
| Zero field critical current $`I_o(mA)`$ | 33 | 2.5 |
| Maximum critical current $`I_{\mathrm{max}}(mA)`$ | 39 | 2.7 |
| Gap quasiparticle current step $`\mathrm{\Delta }I_g(mA)`$ | 88 | 5.2 |
| $`I_{\mathrm{max}}/\mathrm{\Delta }I_g`$ | 0.45 | 0.52 |
| Critical current density $`J_c(A/cm^2)`$ | 3050 | 180 |
| Josephson length $`\lambda _J(\mu m)`$ | 6.9 | 28 |
| Normalized mean circumference $`C/\lambda _J`$ | 72 | 18 |
| Quality factor $`V_m(mV)`$ | 49 | 63 |
| Normal resistance $`R_N(m\mathrm{\Omega })`$ | 36 | 610 |
| ZFS1 asymptotic voltage $`(\mu V)`$ | 51 | 53 |

Table I

## VI The experimental setup

In order to vary the quenching time in the broadest possible range, we have realized the experimental setup schematically shown in Fig.2.

A massive $`Cu`$ block held to the sample holder by two thin $`Cu`$ arms was used to increase the system thermal capacity. The chip was mounted on one side of this block and thermally loosely coupled to it by means of a $`1mm`$ thick teflon sheet. On the other side of the $`Cu`$ block, a thermoblock consisting of a $`50\mathrm{\Omega }`$ carbon resistor and two thermometers in order to measure and to, if necessary, stabilize the $`Cu`$ block temperature, was mounted in good thermal contact. Finally a small sized $`100\mathrm{\Omega }`$ resistor, more precisely a surface mount resistor (SMR), was kept in good thermal contact with the chip by means of a small amount of vacuum grease.

This system, due to the two heating elements placed in tight and loose thermal contact with the chip and with the proper choice of the thermal loads, allowed us to perform the sample quenching over two quite different time scales. In fact, by means of the resistor in the thermoblock, a long time scale was achieved by heating the chip through the $`Cu`$ block and the teflon sheet; on the contrary, a short current pulse through the surface mount resistor on the chip, attained much short thermal cycles.

These two completely different quenching techniques provide timescale ranges that do not overlap, leaving a gap between $`0.2`$ and $`1s`$, that would require a third quenching technique to be filled. We stress that just using a single sample holder with smaller heat capacity would not give us access to the same time-scales. We could not get quench times as short as those from the SMR and, even with a very small exchange gas pressure, could not get times as long as those with the larger sample holder on using the mechanical pumps available to us.

The whole system was kept in a vacuum tight can immersed in the $`LHe`$ bath at $`He`$ gas pipeline pressure. The pressure of $`He`$ gas inside the can could be varied in order to modify the heat exchange between the chip and the environment and, in turn, the speed of the sample cooling. A solenoid was wound around the can to provide a strong vertical magnetic field and Helmholtz coils were instead placed inside the can to generate a weak horizontal magnetic field in the barrier plane in order to tune the critical currents of the annular junctions to their maximum values.

The temperature dependence of the junction gap voltage was exploited to monitor the temperature of the junction itself during the thermal cycle. Fig.3a and Fig.3b show the digitally measured $`V_g(t)`$ for sample $`A`$ current biased on the quasiparticle curve at $`17.7mA`$, (i.e., at about one fifth of $`\mathrm{\Delta }I_g`$ at $`T=4.2K`$), for a slow and fast thermal cycle, respectively. In the case of Fig.3a, a $`100mA`$ current was fed to the thermoblock heater for about $`4s`$ in order to increase the junction temperature up to its critical value where the sample IVC becomes a straight line with a slope corresponding to the junction normal resistance $`R_N`$. In the case of Fig.3b a $`20V`$ high and $`4ms`$ wide, voltage pulse was applied at the SMR. In the Figs.3 the time origin is arbitrarily set at the instant we began to feed the heating elements. It is important to observe that the time scale in Fig.3a is about 50 times larger than that of Fig.3b, although the curve shapes are quite similar.

In our samples the current jump at the gap voltage is very steep and, at $`T=4.2K`$, the voltage changes by less than $`1\%`$ when the current is changed from $`10`$ to $`90\%`$ of the total current jump $`\mathrm{\Delta }I_g`$ and by less of $`10\%`$ at $`T=8.5K`$. Furthermore, at $`T=4.2`$and $`8.5K`$ and for this bias current, the junction voltage was equal to $`2.74`$and $`1.0mV`$, respectively. Therefore, assuming that the electrode gap energies are equal $`\mathrm{\Delta }_1(T)=\mathrm{\Delta }_2(T)=\mathrm{\Delta }(T)`$, and that in the $`4.28.5K`$ range we can neglect the thermal gap smearing, the analytical expression found by Thouless for the gap energy in a strong-coupling superconductor

$$\frac{\mathrm{\Delta }(T)}{\mathrm{\Delta }(0)}=\mathrm{tanh}\frac{\mathrm{\Delta }(T)}{\mathrm{\Delta }(0)}\frac{T_c}{T},$$
(5)
also applies to the junction gap voltage that is proportional to it. The dashed lines in Figs.3 indicate the voltage threshold above which Eq.5 can be used to relate the junction voltage at its temperature, without any significant loss of accuracy. We like to mention that an experimental proof of Eq.5 in $`Nb/Nb`$ tunnel junctions was first evidenced by Broom.

Fig.4a and Fig.4b show the data reported in Fig.3a and Fig.3b, respectively, transformed according to Eq.5, assuming for $`\mathrm{\Delta }(0)`$ and $`T_c`$ the values $`2.85meV`$ and $`8.95K`$, respectively, as found by Monaco et al.  on similar JTJs. Now the dashed horizontal lines indicate the temperature threshold below which the temperature time dependence can be reliably accounted for by our measured data. The large noise at low temperature is the result of an amplification effect of Eq.5, according to which the temperature variation corresponding to a given energy gap variation becomes larger and larger as the temperature becomes smaller. It would be very complicated to write the proper boundary conditions for the heat diffusion equation that would correctly model the full time dependence of the junction temperature. However, for our purposes we are only interested to the cooling process, and we successfully fit our data by a simple thermal relaxation equation:

$$T(t)=T_{fin}+(T_{in}T_{fin})\mathrm{exp}\left(\frac{tt_0}{\tau }\right)$$
(6)
with only two fitting parameters $`t_0`$ and $`\tau `$, being $`T_{in}`$ and $`T_{fin}`$ fixed at $`8.95`$ and $`4.15K`$, respectively. In Eq.6 $`t_0`$ is the time at which $`T=T_{in}=T_c`$ and $`\tau `$ is the relaxation time which sets the cooling time scale. The fitting curves are shown by the thick dashed lines in Fig.4a and Fig.4b, and correspond to a thermal relaxation time $`\tau `$ equal to $`3.6`$ and $`0.073s`$, respectively. The quenching time $`\tau _Q`$ can be obtained from its definition:

$$\frac{T_C}{\tau _Q}=\frac{dT}{dt}_{T=T_C}$$
(7)
giving $`\tau _Q=\tau T_C/(T_{in}T_{fin})`$. For a thermal relaxation from the junction critical temperature down to the helium bath temperature we get $`\tau _Q1.9\tau `$. Eventually, Fig.5 display the values of the quenching times obtained with the process described above, both for the fast (black squares with right vertical scale) and slow (red circles with left vertical scale) cooling processes and for different values of the $`He`$ pressure inside the can. We observe that, by changing the exchange gas pressure and using the two techniques, the quenching time can be changed over a quite large range starting from tenths to tens of seconds. At the end of this paper we will discuss how it is possible to extend this range in both directions; however, as we will show in next section, this range has shown to be large enough for our purposes.

## VII The measurements

Quenching experiments were carried out in a double $`\mu `$-metal shielded cryostat and the transitions from the normal to the superconducting states were performed with no current flowing in the heaters and the thermometers. Both the junction voltage and current leads were shorted during all the thermal cycle. Furthermore, the heat supplied to the sample was such that the maximum temperature reached by the junction was made slightly larger than its critical temperature, say at about $`10K`$, in order to make sure that also the bulk electrode critical temperature ($`T_C9.2K`$) was overcome. In this case, according to Eqs.6 and 7, the value of the quenching time results in a correspondingly smaller $`\tau _Q1.7\tau `$. Due to the approximation made by using Eq.6 and to the experimental uncertainty in the knowledge of the maximum temperature during each thermal cycle, the value of the quenching times has been determined with an overall accuracy as large as $`5\%`$. For each value of the quenching time, in order to estimate the trapping probability, we have carried out a set of 300 thermal cycles and at the end of each cycle the junction IVC was inspected in order to ascertain the possible spontaneous trapping of one or more fluxons.

As we shall see later, the AJTJs are such that the ZK causal length $`\overline{\xi }>C`$ by an order of magnitude when $`\tau _Q=1s`$. Increasing and decreasing $`\tau _Q`$ by an order of magnitude changes $`\overline{\xi }`$ by less than a factor of two. Thus the probability of finding a single fluxon after a quench is small. In the following we will focus our attention only on the probability $`P_1`$ to trap just one fluxon, although a few times we found clear evidence of two and, more seldom, three homopolar fluxons spontaneously trapped during the N-S transition. However, these events were too rare to be statistically significant.

Experimentally, we define $`P_1`$ as the ratio between the number of times in which at the end of the thermal cycle the junction IVC looks like that shown in Fig.1b, i.e. with a tiny critical current and a large first ZFS, and the number of attempts. It is worth to mention here that in the case of simultaneous trapping of a fluxon and an antifluxon, they would annihilate and leave no track of their formation. Therefore, our definition of $`P_1`$ is not rigorous, but it is reasonable as far as the chances to trap two (homo or heteropolar) fluxons are negligibly small. Similarly, an IVC similar to that shown in Fig.1b could be the result of the simultaneous trapping of two fluxons and one antifluxon, or the other way around, but this event is less likely in our experimental situation. For the sake of completeness, it must be added that, in some cases, the IVC displayed either a depressed critical current without ZFS structures or a ZFS with an enhanced critical current. We explain them as due to the trapping of Abrikosov vortices in the junction electrodes and nearby the barrier and we did not take in any account the occurrence of such events, since it is not known if and how the vicinity of Abrikosov vortices influences the fluxon formation.

## VIII The results

When $`\overline{\xi }>C`$ we estimate the probability of finding a fluxon in a single quench to be

$$P_1\frac{C}{\overline{\xi }}=\frac{C}{\xi _0}\left(\frac{\tau _Q}{\tau _0}\right)^\sigma ,$$
(8)
where, from Eq.(4), $`\sigma =0.25`$.

Fig.6 shows on a log-log plot the measured probability $`P_1`$ of a single fluxon trapping obtained by quenching the sample $`A`$ 300 times for each value of the quenching time $`\tau _Q`$ changed by varying the exchange gas pressure and by using both the fast and slow quenching techniques. We observe that the points are quite scattered meaning that the data are statistically poor. Further, for the reasons given earlier there is a gap between fast and slow quenches. Nonetheless, we have clear evidence that i) the trapping of a fluxon occurs on a purely statistical basis being identical the conditions for each thermal cycle and ii) the probability to trap one fluxon is larger when the transition is performed at a faster speed (smaller quenching time) in a qualitative accordance with the causality principle. More precisely, we can distinguish the point to the lower right of the graph corresponding to the slow cycle process which gives in the average a probability to spontaneously trap a fluxon once every fourteen attempts, and the cloud of data to the upper left corresponding to the impulsive junction heating and giving an average probability of one successful event every about six attempts. This suggests that possible temperature gradients induced by the SMR are not an important source of systematic error, since such gradients, with their slowly moving profiles, have a tendency to reduce defect production.

Regardless of the data spread, as suggested if Eq.8 holds true, we attempted to fit the data with an allometric function $`P_1=a\tau _Q^b`$ with $`a`$ and $`b`$ being free fitting parameters. We found that the best fitting curve, shown by the solid line in Fig.6, has a slope $`b=0.27\pm 0.05`$. Such a value of $`b`$, although affected by a $`20\%`$ uncertainty, is in good agreement with the forth root square dependence expected for symmetric junction.

For the coefficient $`a`$ we found the best fitting value of $`0.1\pm 10\%`$ ($`\tau _Q`$ in seconds). This is to be compared with the predicted value of $`C\tau _0^{1/4}/\xi _0`$. Sample $`A`$ had a circumference $`C=500\mu m`$. Its effective superconductor thickness was $`d_s250`$ $`nm`$. At the final temperature $`T_{fin}=4.2K`$, the critical current density was $`J_c(T_{fin})=3050A/cm^2`$ and the Josephson length was $`\lambda _J(T_{fin})=6.9\mu m`$. From this, and $`c_0`$ given earlier, we infer that $`\xi _03.8\mu m`$ and $`\tau _00.17ps`$. This then gives $`C\tau _0^{1/4}/\xi _00.08s^{1/4}`$, in good agreement with the experimental value of $`b`$, given the fact that we only expect agreement in overall normalization to somewhat better than an order of magnitude level. At this level, such a result is immune to systematic error in one or other of the measurement processes. After the problems (discussed below) of the experiments discussed in  to find (reliable) defects at expected densities, if at all, our experiment shows that the ZK estimate remains sensible.

Similar measurements have been carried out for sample $`B`$. Although not in contradiction with (1), the results were affected by a data scattering even larger than that found for sample $`A`$ (shown in Fig.6). This is due to a much smaller normalized length which, according to (8), translates in a expected probability $`P_{1\text{ }}`$, for a given $`\tau _Q`$, about four times smaller (since $`C\tau _0^{1/4}/\xi _00.02s^{1/4}`$ for this sample), far too small too get statistically significant data in reasonable times considering that these measurements are both very time and $`LHe`$ consuming. In order to have data comparable with those of sample $`A`$, sample $`B`$ would have required a quenching times $`4^4`$ times larger. However the roughly measured probability $`P_{1\text{ }}`$of 1 fluxon every 50-100 attempts is in fairly good agreement with the expected value. This shows that our fluxons are not spurious byproducts of the measurement mechanism.

## IX Comments, Future Experiments and Conclusions

We consider this experiment to give a strong confirmation of the Zurek-Kibble predictions. We said in the introductory section of this paper that condensed matter experiments to test the ZK predictions had given mixed results, and it is interesting to put this experiment in that context.

Prior to our experiment, seven other experiments had been performed to test (1), five with fixed $`\tau _Q`$, two with variable $`\tau _Q`$. \[ In addition, the experiment cited earlier on JTJs by two of us (R.M and J.M) was compatible with (4), although it had not been performed with a test of (4) in mind. It was this that motivated the experiment described here.\]

Of those experiments with fixed $`\tau _Q`$, two were experiments on superfluid $`{}_{}{}^{3}HeB`$, which rely on the fact that, when it is bombarded with slow neutrons, energy is released which leads to a hot spot, with temperature $`T>T_c`$, in the superfluid which then cools below $`T_c`$. This leaves behind a tangle of vortices, the topological defects in this system, whose density can be measured. Since $`\tau _Q`$ is fixed by the nuclear process, it is not possible to confirm the predicted value $`\sigma =1/4`$. However, with only a single data point conflating both normalisation and $`\sigma `$ both experiments are highly compatible with (1).

The remaining experiments with fixed $`\tau _Q`$ were two on superfluid $`{}_{}{}^{4}He`$, and one on high temperature superconductors (HTSC).

In principle, the $`{}_{}{}^{4}He`$ experiments, which use a pressure quench with a varying timescale $`\tau _Q`$ to implement the transition, could have allowed for a more complete test, in this case to confirm $`\sigma =1/3`$ (after renormalisation group rescaling). Yet again, vortices are the relevant defects. In practice, the most reliable published experiment sees no vortices. In this context, the vortices seen in an earlier $`{}_{}{}^{4}He`$ experiment, at levels compatible with (1), were most likely an artefact of the experimental setup. Further experiments on $`{}_{}{}^{4}He`$ are underway.

The fifth experiment, on HTSC, measures total flux through a surface i.e. the variance in the topological charge, carried in this case by the Abrikosov vortices. The vortex separation of (1) can be converted into a prediction for the variance, but no flux is seen in contradiction with this prediction, despite the phase separation that is a prerequisite for the result being seen elsewhere. There is no obvious explanation of this null result. An attempt to take gauge fields into account fully shows that there is an additional mechanism for vortex production in the thermal fluctuations of the magnetic field but, as yet, this seems insufficient to explain the result. Such a mechanism will not apply to the JTJs considered in our experiment.

These early experiments have either provided one data point for (1 ), or have been null for whatever reason. Two subsequent experiments have permitted varying quench rates and so an estimate for $`\sigma `$. The most recent involves the B$`\stackrel{´}{e}`$nard-Marangoni conduction-convection transition, in which a homogeneous conduction state is broken into an hexagonal array of convection lines on heating. The defects here are not associated with the line zeroes of an order parameter field, and the viscosity-dependent $`\sigma `$ does not match the ZK prediction, most likely for that reason. The more relevant experiment is carried out in a non-linear optical system, with complex beam-phase the order parameter, satisfying a time-dependent Ginzburg-Landau equation with drift. There has been much numerical analysis of time-dependent Ginzburg-Landau systems, which show agreement with the ZK predictions for scaling exponents. The control parameter in  is not the temperature, but the light intensity. Increasing it leads to pattern formation (defects) at a critical value. The predicted scaling parameter $`\sigma =1/4`$ is recovered to good accuracy as $`\sigma _{exp}=0.25\pm 0.02`$, but agreement with normalisation is not stated.

Given this relatively poor success rate in confirming (1) we are considering a further experiment to measure the ZK scaling exponent, this time with manifestly non-symmetric AJTJs. In  and  we observed that it should have been difficult to make JTJs truly symmetric, as those used here. However, in  and  we had not appreciated how the difference between symmetric and marginally non-symmetric JTJs is smeared by the proximity effect of the $`Al`$ within the insulating layer. Significantly non-symmetric JTJs require different fabrication techniques, but the value of $`\sigma `$ inferred from the same causal arguments is $`\sigma =1/7`$, very different from the value of $`1/4`$ that we tested above. The data from our experiment is incompatible with $`\sigma =1/7`$. This does suggest that a further experiment, with markedly non-symmetric JTJs, should be performed.

Our experiments have demonstrated that quenching time of the order of 1 second give a rather large probability to trap one fluxons on AJTJs having a very large normalized length. However, very long junctions mean very large critical current densities that, in turn, require Josephson barriers so thin that their quality and uniformity is often spoiled; furthermore, in most cases, applications require intermediate length junction or even small junction. For these reasons, it would be highly desirable to compensate the reduced junction length with an increased quenching rate, as it is suggested by the findings for sample $`B`$. Therefore, we like to conclude the paper with some comments on the possible technical improvements that would allow to test the ZK predictions over a broader quenching time range. Firstly, $`\tau _Q`$ can be trivially heightened by increasing the $`Cu`$ block thermal capacitance. On the contrary in order to lower the quenching time, that is to make the $`NS`$ transition faster, it is needed to resort to new techniques since the maximum power that can be dissipated by the surface mount resistors sets an obvious lower threshold on $`\tau _Q`$. One possible way to reach this goal is to perform the junction thermal cycle by means of light pulses. Light dissipates inside the superconducting electrodes, but not in the substrate providing a local junction heating that will relax much faster to the background temperature. We estimate that, by using a properly focussed pulsed light beam, the quenching time scale can be reduced to the $`\mu s`$ range.

## ACKNOWLEDGMENTS

The authors thank L. Filippenko for the sample fabrication and V.P. Koshelets for useful discussions. R.R. thanks the University of Salerno for hospitality. This work is also supported by the COSLAB programme of the European Science Foundation.

## X Table and Figure Captions

* Geometrical and electrical parameters of two selected annular Josephson tunnel junctions at 4.2 K.
* Low voltage part of the experimentally measured current-voltage characteristics of the same annular Josephson tunnel junction a) without trapped fluxons, b) with one trapped fluxon, and c) with two trapped fluxons. For each current branch the corresponding number of travelling fluxons $`F`$ and fluxon-antifluxon pairs $`F\overline{F}`$ is indicated.
* Sketch (dimensions are not to scale) of the cryogenic insert developed to perform the junction thermal cycles with a time scale changing over a broad range.
* Digitally measured time dependence of the junction gap voltage during the a) ”slow” and b) ”fast” thermal cycle. For these measurements the junction was biased at about 1/5 of the gap quasiparticle current step. The horizontal dashed lines indicate the voltage threshold above which Eq.5 can be used to relate the junction voltage at its temperature.
* Time dependence of the junction temperature during the a) ”slow” and b) ”fast” thermal cycle obtained from the data of Figs.3 transformed according to Eq.5. The horizontal dashed lines indicate the temperature threshold below which the temperature time dependence can be reliably accounted for by our measured data. Furthermore, the thick dashed lines are best fitting curve of the cooling process according to the thermal relaxation expression Eq.6.
* Quenching time $`\tau _Q`$ as a function of the $`He`$ gas pressure inside the vacuum tight can. By changing the pressure of the exchange gas inside the can, the system thermal constants are varied. The solid squares refer to the left vertical scale while the solid circles refers to the right vertical scale.
* Log-log plot of the measured probability $`P_1`$ to trap one fluxon versus the quenching time $`\tau _Q`$. The solid line is the best fitting curve found assuming a power law dependence as suggested by Eq.1. To a good degree of approximation the fit is in agreement with a forth root square dependence as expected for symmetric annular Josephson tunnel junctions.


text:
# Sulphur and iron abundances in halo stars1

## 1. Introduction

Sulphur is generally regarded as an $`\alpha `$-capture element. The work on Galactic stars by François (1987, 1988) supported this view by showing that \[S/Fe\] increases from zero at solar metallicities to a plateau level of about +0.5 dex in the metallicity range $`1.8<\text{[Fe/H]}<0.8`$, an analogous behaviour to that of other $`\alpha `$-elements Mg, Si, and Ca, see Ryan et al. (1996). The standard interpretation is that this trend arises from the time delay in the production of 2/3 of the iron by SN of Type Ia relative to the near-instantaneous release of the $`\alpha `$-elements by Type II SN.

Recent observations of sulphur in metal-poor stars by Israelian & Rebolo (2001) have, however, challenged this view. Their data suggest that \[S/Fe\] increases linearly with decreasing \[Fe/H\] to a level as high as $`\text{[S/Fe]}+0.7`$ at $`\text{[Fe/H]}=2.0`$. The study of Takada-Hidai et al. (2002) based on Keck HIRES observations supports a quasi linear dependence of \[S/Fe\] on \[Fe/H\] although in their case \[S/Fe\] reaches only 0.5 dex at $`\text{[Fe/H]}=2.0`$.

As a possible explanation of the high value of \[S/Fe\] in metal-poor stars, Israelian & Rebolo propose that massive supernovae with exploding He-cores and a high explosion energy make a significant contribution to the early chemical evolution of galaxies. According to Nakamura et al. (2001) these hypernovae overproduce S with respect to O, Mg and Fe. With this intriguing possibility in mind a more thorough investigation of sulphur abundances in halo stars seems worthwhile. A clarification of the trend of S abundances is also much needed in deciphering the chemical enrichment of damped Ly$`\alpha `$ systems (DLAs), widely regarded as the progenitors of present-day galaxies at high redshift. Its importance stems from the fact that, unlike most other heavy elements, S is not depleted onto dust. Consequently, observations of the relatively weak S ii triplet resonance lines at $`\lambda \lambda 1250,1253,1259`$ yield a direct measurement of the abundance of S in DLAs. The only other element for which this is the case is Zn, and indeed most of our current knowledge of the chemical evolution of the universe at high redshift is based on surveys of \[Zn/H\] in DLAs (e.g. Pettini et al. 1999, Prochaska & Wolfe 2002). If S is an $`\alpha `$-capture element, then its abundance relative to Zn (an iron-peak element) could be used as ‘a chemical clock’ to date the star-formation process at high $`z`$. Specifically, if the major star formation episodes in the DLAs observed occurred within the last $`0.5`$ Gyr, we would expect to measure enhanced \[S/Zn\] ratios, and vice versa. Data on the abundance of S in DLAs have been relatively scarce until recently, but are now becoming available at a progressively faster rate thanks largely to the high sensitivity of the VLT/UVES high resolution spectrograph at blue and near-UV wavelengths. The picture is still confused (Centurión et al. 2000, Prochaska & Wolfe 2002), but one thing is clear: without a secure knowledge of the behaviour of S in metal-poor Galactic stars we stand no hope of interpreting the situation at high $`z`$.

In this paper we report on a large scale survey of sulphur abundances in metal-poor stars carried out with VLT/UVES. Effects on the derived S/Fe ratio from the modelling of stellar atmospheres are discussed, and preliminary results for about half of the 35 stars observed are presented.

## 2. Observations and data reduction

Previous determinations of sulphur abundances in metal-poor stars, including those of Israelian & Rebolo (2001) and Takada-Hidai et al. (2002), have been based on the high excitation ($`\chi =7.87`$ eV) $`\lambda \lambda `$8693.9, 8694.6 S i lines. They are, however, very weak in metal-poor dwarfs and giants ($`W<`$ 3 mÅ at $`\text{[Fe/H]}2.0`$), and practically impossible to observe at $`\text{[Fe/H]}2.5`$. In our survey we have taken advantage of the high efficiency of the VLT/UVES instrument and concentrated on the stronger S i triplet ($`\chi =6.52`$ eV) at 9212 - 9238 Å. The dichroic mode of UVES was used to cover the spectral region 6700 - 10000 Å in the red arm and 4000 - 5000 Å in the blue arm; the latter includes a number of Fe ii lines suitable for determining the iron abundance. The resolution of the spectra is $`\lambda /\mathrm{\Delta }\lambda \mathrm{60\hspace{0.17em}000}`$ with 4 pixels per spectral resolution element. Typical $`S/N`$ ratios are 300 in the blue region and 200 at the 9212 - 9238 Å triplet.

The red spectra were reduced using the echelle package in IRAF. A problem with the 9212 - 9238 Å region is the presence of numerous, quite strong telluric H2O lines, see Fig. 1. To remove these lines fast rotating, early-type stars were observed on each night and reduced in the same way as the program stars. The IRAF task ‘$`\mathrm{t}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{u}\mathrm{r}\mathrm{i}\mathrm{c}`$’ was then used to remove the telluric lines. When the early-type star is observed at about the same airmass as the program star this technique also serves to remove a residual fringing remaining from the flatfielding of the spectra.

Representative spectra of stars are shown in Fig. 2 (8694 Å region) and Fig. 3 (9212 Å region). As seen the S i line at 9212.9 Å can be clearly detected in a $`\text{[Fe/H]}3.2`$ star. Eight stars were observed twice on different nights; a comparison shows that the equivalent widths of the $`\lambda \lambda `$8693.9, 8694.6 lines are measured with an accuracy of about $`\pm 0.5`$ mÅ whereas the $`\lambda \lambda `$9212.9, 9237.5 lines are measured to an accuracy of $`\pm 2`$ mÅ.

The blue spectra were reduced with MIDAS routines. For a few stars it was checked that IRAF reduction gives practically the same result. These blue spectra have $`S/N300`$ and allowed us to measure accurate equivalent widths of about 20 weak Fe ii lines.

## 3. Modelling of the stellar atmospheres

The measured equivalent widths were used to derive S and Fe abundances based on 1D MARCS (Asplund et al. 1997) model atmospheres assuming LTE. Effective temperatures were derived from the color indices $`by`$ and $`VK`$, and surface gravities via Hipparcos parallaxes and/or the Strömgren $`c_1`$ index. For details we refer to the recent paper by Nissen et al. (2002) on oxygen abundances.

Oscillator strengths of the blue Fe ii lines were determined from an inverted abundance analysis of three stars (HD 103723, HD 160617 and HD 140283) adopting \[Fe/H\] values determined by Nissen et al. (2002) from very accurate equivalent widths of weak Fe ii lines analyzed differentially with respect to the solar spectrum. This ensures that we are on the same \[Fe/H\] scale as Nissen et al. (2002).

Oscillator strengths of the sulphur lines were taken from Lambert & Luck (1978). We note that the values for the 8693.9 and 8694.6 Å lines (log$`gf`$ = $`0.56`$ and 0.03, respectively) and the solar spectrum equivalent widths provide a sulphur abundance of the Sun very close to the meteoritic value, log$`ϵ`$(S) = 7.20. The 9212.9 and 9237.5 Å S i are too strong in the solar spectrum to provide reliable abundances, but we note that the sulphur abundance derived from this pair of lines agree well with that derived from the $`\lambda \lambda `$8693.9, 8694.6 pair for 12 stars in the metallicity range $`1.5<\text{[Fe/H]}<0.7`$. The mean difference between the S abundance derived from the $`\lambda \lambda `$8693.9, 8694.6 pair and the 9212.9,9237.5 Å lines is +0.03 dex with an rms deviation of the difference of $`\pm 0.10`$ dex. The third line of the near-infrared S i triplet at 9228.1 Å is unfortunately very close to the Paschen $`\zeta `$ H i line at 9229.0 Å and could not be used for a reliable determination of the sulphur abundance.

Errors in $`T_{\mathrm{eff}}`$ ($`\pm 70`$ K) and log $`g`$ ($`\pm 0.15`$ dex) have only a minor effect on the derived \[S/Fe\] (less than 0.05 dex). This is connected to the fact that the abundances of both elements are determined from their main ionization stages, S i and Fe ii, respectively. The microturbulence was determined from requiring that the abundance derived from the 20 Fe ii lines should be independent of $`W`$, or in the case of stars with very weak Fe ii lines ($`\text{[Fe/H]}<1.5`$) $`\xi _{\mathrm{micro}}`$ is assumed to be 1.5 km s-1. In any case, the effect of the uncertainty of $`\xi _{\mathrm{micro}}`$ on the S abundance is quite negligible due to the weakness of the S i lines.

As shown by Takada-Hidai et al. (2002) non-LTE effects on the weak, high excitation S i lines are likely to be small ($`<0.05`$ dex). The same is the case for the Fe ii lines according to Thévenin & Idiart (1999).

Much attention has recently been paid to the effect on abundance determinations of 3D hydrodynamical modelling of convection and granulation in stellar atmospheres. As first shown by Asplund et al. (1999) there is a very significant surface cooling in metal-poor stars caused by the adiabatic cooling of rising and expanding elements coupled with the lack of radiative heating due to the low line opacity in metal-poor stars. The effect on lines formed in the upper atmosphere can be very large. As an example Asplund & García Pérez (2001) estimated 3D effects on oxygen abundances derived from OH lines of the order of $`0.6`$ dex for turnoff stars with $`\text{[Fe/H]}<2.0`$.

To investigate the 3D effects on the sulphur lines used in the present study a differential 1D - 3D study similar to that described by Asplund & García Pérez (2001) has been carried out. The results are given in Table 1 for two sets of model atmospheres with metallicities of $`\text{[Fe/H]}=0.0,1.0,2.0`$ and $`3.0`$. The first set has $`T_{\mathrm{eff}}`$ and log $`g`$ close to the solar values, whereas the second set represents stars close to the turnoff. For each model atmosphere Table 1 first gives the adopted S abundance, then the corresponding computed 1D equivalent width, and finally the S abundance calculated from the 1D equivalent width on the basis of the 3D model atmosphere. As seen the maximum 3D effect occurs for the most metal-poor stars: the 3D abundance is increased by 0.07 to 0.13 dex with respect to the 1D abundance in the case of the $`\lambda \lambda `$9212.9, 9237.5 S i lines. Interestingly, there is, however, a similar 3D effect ($`+0.10`$ dex) on the Fe abundance derived from Fe ii lines (Nissen et al. 2002). Hence, it appears that we are in the favourable situation that the S/Fe ratio is quite immune to 1D - 3D effects.

## 4. Results

Up to now we have reduced and analyzed about half of the VLT/UVES spectra. The derived sulphur abundances for 19 halo stars are given in Table 2 together with atmospheric parameters and iron abundances. Fig. 4 shows a plot of \[S/Fe\] vs. \[Fe/H\] with disk stars from Chen et al. (2002) added. The sulphur abundances of the disk stars were derived from the $`\lambda \lambda `$8693.9, 8694.6 S i lines and three other weak S i lines in the spectral range 6000 – 6800 Å.

As seen from Fig. 4 \[S/Fe\] shows a clear plateau ($`\text{[S/Fe]}0.35`$) for \[Fe/H\] below $`1.2`$. There is no tendency for a linear rise of \[S/Fe\] as claimed by Israelian & Rebolo (2001) and Takada-Hidai et al. (2002). We note in this connection that their samples include only three stars with $`\text{[Fe/H]}<2.0`$ and none below $`\text{[Fe/H]}=2.3`$, whereas our sample includes seven stars in the range $`3.2<\text{[Fe/H]}<2.0`$. Furthermore, the error bars on their sulphur abundance determinations are large due to the weakness of the S i 8694.6 Å line; the 9212.8 and 9237.5 Å S i lines, on the other hand, are rather ideal for sulphur abundance determinations below $`\text{[Fe/H]}=2.0`$. Although we cannot exclude the existence of a few stars with high values of S/Fe222We do confirm the high value $`\text{[S/Fe]}0.7`$ of HD 19445 ($`\text{[Fe/H]}1.9`$) if we adopt the equivalent width (2.7 mÅ) of the S i 8694.6 Å line given by Israelian & Rebolo based on their very high S/N spectrum it would appear that sulphur shows the same trend as other typical $`\alpha `$-capture elements Mg, Si and Ca. Hence, our preliminary conclusion is that there is no need to invoke element production by hypernovae or very massive supernovae to explain the general behaviour of sulphur. Traditional Galactic evolution models with near-instantaneous production of $`\alpha `$-elements by Type II supernovae and delayed production of the iron-peak elements are in good agreement with our data for sulphur (Chiappini et al. 1999, Goswami & Prantzos 2000).

An interesting detail in Fig. 4 should be noted. All stars from the present investigation (plotted with filled circles) have halo kinematics, including six stars with $`0.9<\text{[Fe/H]}<0.7`$. They have Galactic rotational velocities less than 50 km s-1 i.e. well below the characteristic rotational velocity of disk stars (225 km s-1) and thick disk stars (175 km s-1). As seen from the figure, there is then an overlap between halo stars and thick disk stars in the metallicity range $`1.0<\text{[Fe/H]}<0.6`$. Out of the six stars four have enhanced S/Fe ratios like the thick disk stars but two (HD 103723 and HD 105004) have a solar S/Fe ratio. Interestingly, HD 103723 and HD 105004 were also found to be deficient in other $`\alpha `$-element (O, Mg, Si, Ca, Ti) by Nissen & Schuster (1997). On the basis of their Galactic orbits, these authors speculated that the $`\alpha `$-deficient stars have been accreted from dwarf galaxies with a slow chemical evolution in contrast to the ‘normal’ halo stars with enhanced \[$`\alpha `$/Fe\] formed in a fast chemical evolution in the inner part of our Galaxy. The fact that S and the classical $`\alpha `$-elements, Mg, Si and Ca, are deficient in the same stars is a further indication that sulphur belongs to this group of elements. Thus there seems to be little ground, on the basis of our data, for the reservations expressed by Prochaska et al. (2000) concerning use of S as an $`\alpha `$-element in the analysis of abundance ratios in DLAs. On the other hand, one should remain aware that the overabundance of the $`\alpha `$-elements at low metallicities need not be a universal phenomenon. The finding by Nissen & Schuster (1997) that even some stars in the Milky Way do not fit this pattern is a warning against the naive assumption that such an overabundance should also apply to all metal-poor DLAs.

## References

Asplund M., Gustafsson B., Kiselman D., Eriksson K. 1997, A&A, 318, 521

Asplund M., Nordlund Å., Trampedach R., Stein R.F. 1999, A&A, 346, L17

Asplund M., García Pérez A.E. 2001, A&A, 372, 601

Centurión M., Bonifacio P., Molaro P., Vladilo G. 2000, ApJ, 536, 540

Chen Y.Q., Nissen P.E., Zhao G., Asplund M. 2002, A&A in press, (astro-ph/0206075)

Chiappini C., Matteucci F., Beers T.C., Nomoto K. 1999, ApJ, 515, 226

François P. 1987, A&A, 176, 294

François P. 1988, A&A, 195, 226

Goswami A., Prantzos N. 2000, A&A, 359, 191

Israelian G., Rebolo R. 2001, ApJ, 557, L43

Lambert D.L., Luck R.E. 1978, MNRAS, 183, 79

Nissen P.E., Schuster W.J. 1997, A&A, 326, 751

Nissen P.E., Primas F., Asplund M., Lambert D.L. 2002, A&A in press, (astro-ph/0205372)

Pettini M., Ellison S.L., Steidel C.C., Bowen D.V. 1999, ApJ, 510, 576

Prochaska J. X., Naumov S.O., Carney B.W. et al. 2000, AJ, 120, 2513

Prochaska J. X., Wolfe A.M. 2002, ApJ, 566, 68

Ryan S.G., Norris J.E., Beers T.C. 1996, ApJ, 471, 254

Takada-Hidai M., Takeda Y., Sato S., et al. 2002, ApJ, 573, 584

Thévenin F., Idiart T.P. 1999, ApJ, 521, 753

Nakamura T., Umeda H., Iwamoto K.I. et al. 2001, ApJ, 555, 880


text:
# Metastable Random Field Ising model with exchange enhancement: a simple model for Exchange Bias

## I Introduction

Hysteresis and metastability are intriguing phenomena with implications in both fundamental and applied physics. Magnetic systems are the prototypical example of thermodynamic systems exhibiting hysteresis cycles for which different theoretical approaches have been proposed G.Bertotti (1998). Besides the classical micro-magnetic analysis, based on a continuous description of the magnetic properties of the system, more recently a lot of efforts have been devoted to the study of lattice models. For example, the zero-temperature Random Field Ising model (RFIM) driven by an external field with convenient metastable dynamics has been very successful in qualitatively explaining some basic properties of rate independent hysteresis loops Sethna et al. (1993); Vives and Planes (1994); E.Obradó et al. (1999); E.Vives and A.Planes (2001). The most important achievement of this model has been to give a simultaneous explanation of the effect of disorder on the hysteresis loops and the existence of Barkhaussen noise with critical properties. Less attention has been paid to the use of such models for understanding other interesting features of the hysteresis loops such as remanence, coercivity, minor loop properties or exchange bias (EB) J.Nogués and I.K.Schuller (1999); A.E.Berkowitz and K.Takano (1999), which is the property we will focus our attention on here.

We present a modification of the zero-temperature RFIM that allows magnetic hysteresis loops with EB to be reproduced. The main characteristic of EB is that the hysteresis loops, represented as magnetization $`m`$ versus external applied field $`H`$, are not centred on $`H=0`$ but exhibit a displacement in the field axis by an amount $`H_{EB}`$ (exchange bias field) J.Nogués and I.K.Schuller (1999); A.E.Berkowitz and K.Takano (1999). This property has received a lot of attention recently, since the possibility of finding systems with large EB has enormous technological interest B.Dieny et al. (1991). Experimentally EB has been found in different magnetic systems J.Nogués and I.K.Schuller (1999); A.E.Berkowitz and K.Takano (1999). The basic ingredient for EB is the existence of interfaces between ferromagnetic (FM) and antiferromagnetic (AFM) systems, where coupling can be induced after field cooling from above the Neel temperature of the AFM. This heat treatment freezes some of the magnetic moments at the interface which are supposed to be responsible for the occurrence of EB. The prototype is a FM/AFM bilayer, for instance Co/CoO Miltényi et al. (2000), NiFe/NiMn G.Li et al. (2000), Fe/FeF2 J.Nogués et al. (1996) and Fe/MnF2 C.Leighton et al. (2000). This effect is also observed in granular systems formed by small particles with a ferromagnetic core covered by their native antiferromagnetic oxide S.Gangopadhyay et al. (1993).

Different models have been proposed to understand EB. Although a basic qualitative explanation was given 40 years ago W.H.Meiklejohn (1962) a deep understanding of the phenomenon has not yet been achieved M.Kiwi (2001); R.L.Stamps (2000). Different features remain unclear: the role played by the AFM thickness T.Ambrose and C.L.Chien (1996); Miltényi et al. (2000), the formation of domain walls A.P.Malozemoff (1987); D.Mauri et al. (1987), whether the frozen spins belong to the FM or to the AFM layer M.Kiwi et al. (1999), etc… Especially intriguing is the fact that EB not only occurs, in uncompensated AFM layers which exhibit a net magnetization after being cooled, but also in compensated AFM layers with zero net interface magnetization N.C.Koon (1997).

The aim of the present paper is to introduce a very simple model with a new mechanism for the explanation of EB in totally compensated layers. The model is based on a lattice spin system with metastable dynamics for which some of the exchange interactions show a marked enhancement. In section II the Hamiltonian and the detailed mestastable dynamics are presented. In section III we show the results of the numerical simulations. In section IV we discuss the possible physical origin of the exchange enhancement. In section V we compare with available experimental data and, finally, in section VI we summarize and conclude.

## II Model

The model is intended to reproduce the properties of the ferromagnetic layer only, and regards the AFM part to be totally quenched and thus does not contribute to the net magnetization (compensated AFM layer). Consequently, we consider the 2d RFIM on the square lattice, although a generalization to bulk ferromagnets or thin layers could easily be implemented. Note, however, that the AFM plays an indirect role, as will be discussed later.

The mathematical formulation of the model is very similar to the RFIM on a square lattice with size $`N=L\times L`$. On each lattice site we define a spin variable $`S_i`$ which takes values $`\pm 1`$. The Hamiltonian, in reduced units, reads:

$$=\underset{ij}{\overset{nn}{}}J_{ij}S_iS_j\underset{i}{\overset{N}{}}h_iS_iH\underset{i}{\overset{N}{}}S_i$$
(1)
The first sum is the ferromagnetic exchange contribution that extends over nearest neighbour pairs ($`J_{ij}>0`$). The second sum accounts for the interaction with quenched random fields $`h_i`$, which stands for the disorder present in any ferromagnetic system. The last term is the interaction with the external driving field $`H`$. The random fields $`h_i`$ are independent and distributed according to a Gaussian probability density:

$$p(h_i)=\frac{1}{\sqrt{2\pi }\sigma }e^{\frac{h_i^2}{2\sigma ^2}}$$
(2)
where $`\sigma `$ is the standard deviation of the random fields and controls the amount of disorder in the system. The novelty of the model is in the values of the exchange constants $`J_{ij}`$, which are not equal for all spin pairs: we consider that $`J_{ij}=J`$ except for a fraction $`f`$ of the bonds (selected at random) for which $`J_{ij}=J_E>>J`$. This fraction of bonds is supposed to contain the effect of the quenched antiferromagnetic layer. A physical reason for this local exchange enhancement will be discussed in section IV. Such a bond distribution can be mathematically expressed as:

$$p(J_{ij})=(1f)\delta (J_{ij}J)+f\delta (J_{ij}J_E)$$
(3)
We have focused our study in the region of small values of $`f`$ ($`f<0.06`$). The magnetization of the system is defined as $`m=_{i=1}^NS_i/N`$. For the analysis of the hysteresis loops we use the so called synchronous local relaxation dynamics. This is the standard dynamics used in previous studies of the zero-temperature RFIM Sethna et al. (1993). Each spin $`S_i`$ flips according to the sign of its local field $`H_i`$ given by:

$$H_i=\underset{j=1}{\overset{4}{}}J_{ij}S_j+H+h_i$$
(4)
where the first sum extends over the four neighbours of $`S_i`$. We start with a value of $`H`$ large enough so that the stable situation is given by all the spins $`S_i=1`$. We decrease the external field $`H`$ until $`H_i`$ vanishes on a certain spin. The spin is then reversed keeping $`H`$ constant. This reversal may destabilize some of the neighbouring spins which are then reversed simultaneously. This is the beginning of an avalanche. The avalanches proceed until a new stable situation with all the spins $`S_i`$ aligned with their respective local fields $`H_i`$ is reached. We can then continue decreasing the external field $`H`$.

Most of the calculated properties are averaged over a large number ($`10^3`$) of different realizations of disorder. Averages are indicated by the symbol $``$. We will consider, without loss of generality, that $`J=1`$. Therefore, from now on, magnetic fields and energies are given in units of $`J`$.

## III Results

In Fig. 1 we present an example of a hysteresis loop obtained with the numerical simulation of a system with $`L=50`$, $`J_E=20`$, $`f=0.03`$ and $`\sigma =1.65`$. The external field is cycled between $`H=\pm 2.7`$. As can be seen, the loop exhibits remarkable EB.

At first glance it may look surprising that the model defined in the previous section displays such asymmetry, since the Hamiltonian is totally symmetric under the changes $`S_iS_i`$ and $`HH`$. The reason is that the hysteresis loop shown in figure 1 corresponds, strictly speaking, to a minor loop. This is revealed in Fig. 2 where the system is cycled between $`H=\pm 19.4`$ (which is a field that is one order of magnitude larger than the coercive field $`H_{c1}`$). Due to the existence of a tiny fraction of very large exchange interactions, the total loop exhibits long, flat plateaux in which the system behaves reversibly exactly as if it was saturated. Only when cycling between extremely large values of the external field, does one obtain the symmetric hysteresis loop. Therefore, the loops with EB are incomplete loops and are accompanied by a magnetization shift.

In order to perform a quantitative analysis of this behaviour, averaging over different realizations of disorder is carried out. The hysteresis loops are systematically obtained according to the following protocol: decreasing the field from $`H=+\mathrm{}`$ to $`2H_{c1}<0`$ and increasing the field again to $`H=+\mathrm{}`$, where $`H_{c1}`$ is the coercive field in the decreasing branch. We also compute the pseudo coercive field $`H_{c2}`$ in the increasing branch (see figure 1).

The criteria for choosing the value $`2H_{c1}`$ as a returning point is similar to the criteria used in many experimental cases. One could easily change this limit to $`3H_{c1}`$ or $`4H_{c1}`$ without changing the results provided that $`J_E`$ is large enough. This can be easily understood from the flat tails in the full hysteresis loop shown in figure 2.

The EB field $`H_{EB}`$, the coercivity $`\mathrm{\Delta }H`$ and the magnetization shift $`m_B`$ are defined as:

$$H_{EB}=\frac{H_{c1}+H_{c2}}{2}$$
(5)
$$\mathrm{\Delta }H=H_{c2}H_{c1}$$
(6)
$$m_B=\frac{m_{s1}+m_{s2}}{2}=\frac{1+m_{s2}}{2}$$
(7)
According to this definition, the loops shifted to the left on the $`H`$ axis (as occurs with the loops in the present paper) will have negative exchange bias field. Figures 3 and 4 show the dependence of $`|H_{EB}|`$ and $`\mathrm{\Delta }H`$ on $`J_E`$ for $`\sigma =1`$ and different values of $`f`$ as indicated by the legend. Even for very low values of $`f`$, $`|H_{EB}|`$ increases and saturates for large enough values of $`J_E`$.

In the case of coercivity (Fig. 4) two important results should be emphasized: the increase (almost 40 % in certain cases) in coercivity for intermediate values of $`J_E`$ and the saturation at a constant value (which depends on $`f`$) for large $`J_E`$. Such limiting values at large $`J_E`$ however, are smaller than the coercivity of the system without exchange enhancement.

In order to analyze the dependence of the system properties on the amount of disorder $`\sigma `$, we choose a value of $`J_E`$ that is large enough so that $`|H_{EB}|`$ has reached the limiting maximum value (see Fig. 3). Figures 5 and 6 show the behaviour of $`|H_{EB}|`$ and $`\mathrm{\Delta }H`$ versus $`\sigma `$ for $`J_E=20`$ and different values of $`f`$ as indicated by the legend. Unexpectedly, $`|H_{EB}|`$ shows non-monotonic behaviour with $`\sigma `$, first decreasing until a minimum is reached, but which increases slowly for large amounts of disorder.

The marked variation of $`|H_{EB}|`$ and $`\mathrm{\Delta }H`$ for values of $`f`$ between 0.015 and 0.025 when $`\sigma 0`$ is associated with the fact that the ascending part of the hysteresis loop is very sensitive to the existence of nuclei of unreversed spins in the negative magnetized state. For very low values of $`f`$ we expect that all the nuclei of unreversed spins will be formed by two positive spins joined by an enhanced bond. The negative spins surrounding such a nucleus will flip (in the $`\sigma 0`$ limit) around $`H=2`$ (see Eq. 4). However, for larger values of $`f`$, larger nuclei will exist. For instance, a nucleus formed by three spins joined by two perpendicular bonds, both with exchange enhancement, acts as a nucleating seed which triggers the avalanches towards the positive magnetization phase when $`H=0`$. When $`f`$ is large enough, such that the probability for such nuclei is significantly different from zero, the coercive field for the ascending branch decreases from 2 to 0, thus increasing $`H_{EB}`$ and decreasing $`\mathrm{\Delta }H`$.

The non-monotonic behaviour of the EB when disorder is increased better seen by plotting $`|H_B|/\mathrm{\Delta }H`$, which is a dimensionless quantity and is more relevant from the experimental point of view. Note that for large $`f`$ and $`\sigma `$, we can find strongly biased hysteresis loops for which $`H_{c1}`$ and $`H_{c2}`$ are negative. An example, obtained by sweeping the field between $`H=\pm 4.5`$ is shown in figure 8.

## IV Discussion

One of the reasons for the present work was to obtain hysteresis loops with EB by modifying the zero-temperature RFIM as little as possible. The straightforward, naive idea would be to consider non-symmetric distributions of random fields, for which the average $`h_i0`$ will create the displacement of the loop. From our point of view, this will correspond to the effect of an uncompensated AFM layer. Constraining ourselves to the inclusion of compensated disorder, we have found that any symmetric distribution of random fields cannot give EB. To understand this, suppose that a certain fraction of spins is pinned by very positive and (symmetrically) very negative random fields. Fig. 9(a) shows the schematic hysteresis loop corresponding to such a system. The spins with more positive random fields, which are the last to reverse in the decreasing branch, will be the first to flip in the increasing branch of the full hysteresis loop. Therefore, the full hysteresis loop (symmetric, without EB) will overlap with the minor hysteresis loop, as shown in Fig. 9a. In contrast, when the bond distribution is distorted, as done in the model presented in this work (and which preserves the symmetry properties of the Hamiltonian in Eq. 1), the spins with larger ferromagnetic coupling will be the last to reverse in the decreasing branch and also the last to be reversed in the increasing branch. Thus, minor loops will not coincide with the full loops and will easily exhibit large EB as indicated in Fig. 9b.

The justification of the proposed model would require a physical explanation for the exchange enhancement phenomena. Recently it has been suggested that spin-waves in the FM/AFM interface could be responsible for an enhancement of the exchange coupling H.Suhl and I.K.Schuller (1998). Here we propose a different mechanism based on the existence of quenched disorder in the AFM layer. The most common case is the existence of antiphase domain boundaries as considered in Refs. Miltényi et al., 2000, U.Nowak et al., 2002, J.Keller et al., 2002 and A.Misra et al., 2001. One can assume that the exchange interaction between the magnetic moments in the FM layer has two contributions: the first comming from the direct overlap of the electronic wave functions of the atoms in the FM layer and a second arising from a superexchange interaction through the overlap with the electrons in the AFM layer. The existence of quenched disorder in the AFM layer can modify this second contribution of the exchange interaction giving rise to the exchange enhancement. Since the energies associated with the broken AFM bonds can be higher than the FM exchange energies (for instance due to the existence of strong anisotropy), it is plausible to imagine that the defects in the AFM can influence the FM exchange interactions.

We would like to give a possible mathematical formulation for such a physical mechanism within the framework of lattice models. Let us consider that the interaction between two neighbouring spins $`S_i`$ and $`S_j`$ of the FM layer is given by:

$$E_{ij}=J_0S_iS_jKS_iS_j\sigma _i\sigma _j$$
(8)
where $`\sigma _i`$ and $`\sigma _j`$ are the spin variables describing the magnetic moments in the AFM layer that sit exactly below the $`S_i`$ and $`S_j`$ spins of the FM layer. The constant $`J_0>0`$ accounts for the exchange interactions in the free FM layer and $`K>0`$ accounts for the coupling between the two layers. This coupling term can include, in an effective form, the interactions with many atomic layers in the AFM. Figure 10 shows a schematic representation of the spins of the FM and the AFM layer in two situations: (a) corresponds to the normal situation for an ordered AFM layer, whereas (b) corresponds to the case in which the AFM exhibits an antiphase domain boundary.

In case (a), the energy of the $`i,j`$ pair is $`E_{ij}=(J_0K)S_iS_j`$, whereas in case (b) (or on any antiphase boundary) the energy is given by $`E_{ij}=(J_0+K)S_iS_j`$. The ratio between the two exchange constants is $`J_E=(J_0+K)/(J_0K)`$, which can be much larger than $`1`$ when $`K`$ is close but smaller than $`J_0`$.

## V Comparison with experiment

The idea that exchange bias could be a minor loop effect was already suggested to explain the influence of the AFM spin-flop in EB J.Nogués et al. (2000a). Moreover, it has been shown experimentally that the uncompensated spins which have been suggested to result in exchange bias K.Takano et al. (1997) and which could result in the proposed exchange enhancement, can be reversed at high enough fields T.Ambrose and C.L.Chien (1996) (much higher than those usually used in exchange bias studies). Hence, in this case EB could also be considered as a minor loop effect to a certain extent.

Moreover, the peak observed in the behaviour of the coercivity with $`J_E`$ (see fig. 4) could be correlated with the experimental observation of peaks in the coercive field $`H_c(t_{AFM})`$ (where $`t_{AFM}`$ is the thickness of the AFM layer) close to the critical AFM thickness for the onset of EB or in $`H_c(T)`$ close to the Néel temperature (see figures 10 and 13 in Ref. J.Nogués and I.K.Schuller, 1999). One could consider that as the AFM thickness or the temperature increases, the FM-AFM coupling, which causes the enhancement of JE, decreases and is consequently the origin of the observed $`H_c`$ increase, in agreement with theory. Similarly, the changes in $`H_{EB}`$ and $`\mathrm{\Delta }H`$ with JE could also be correlated with the behaviour experimentally observed in Fe/MnF2 under a cooling field C.Leighton et al. (2000).

Another interesting result of this model is the behaviour of EB with the amount of disorder $`\sigma `$, as seen in Figs. 5 and 7. Experimentally the role played by disorder in EB is not clear. In some cases there is evidence that increasing disorder increases $`H_{EB}`$, whereas in other cases the opposite effect has been found J.Nogués et al. (1999). Our model is able to explain both possibilities. Moreover, it is noteworthy that experimental non-monotonic dependence of $`H_{EB}`$ on disorder (e.g. roughness, irradiation damage or structural disorder) have also been reported Miltényi et al. (2000); C.Leighton et al. (1999); A.Mougin et al. (2000).

Finally, another remarkable result from our model is that the loops exhibiting EB also exhibit a vertical shift in the magnetization axis. This effect has been observed experimentally W.H.Meiklejohn and C.P.Bean (1957); F.Radu et al. (2002); J.Nogués et al. (1999); J.Nogués et al. (2000b); J.Keller et al. (2002). However, since the fraction $`f`$ is very small, the shift can also be very small and in some cases, would be difficult to observe experimentally. Fig. 11 shows the dependence of the magnetization shift $`m_B`$ on $`H_{EB}`$ for different values of $`f`$. The points correspond to different values of $`\sigma `$ ranging from $`0`$ to $`2`$. In view of these results, it would be interesting to measure such displacements in different experimental systems.

## VI Summary and conclusions

We have presented a simple model for the study of Exchange Bias. The model is based on the RFIM driven by an external field with metastable dynamics. The key ingredient is the existence of a tiny fraction of ferromagnetic bonds, which are strongly enhanced. This creates a long reversible plateaux in the decreasing branch. When this pseudo-saturated state is reached, the reversal of the field gives rise to an EB loop. Different properties of these loops have been computed as a function of the fraction of enhanced bonds and the amount of disorder in the system. We have suggested a possible physical mechanism to justify the existence of such an exchange enhancement in a FM layer on a AFM layer with antiphase domain boundaries, based on the existence of superexchange coupling between the two layers.

The main conclusions of the paper are: (i) EB is due to a minor loop effect; (ii) compensated AFM layers can exhibit exchange biased loops with a concomitant magnetization shift; (iii) many experimental phenomena related to exchange bias such as: peaks in the coercive field, magnetization shifts, marked coercivity increase or non-monotonic dependence of $`H_{EB}`$ on disorder can successfully be reproduced with this model.

## VII Acknowledgements

We acknowledge many fruitful discussions with J.Nogués. This work has received financial support from CICyT (Spain), project MAT2001-3251 and CIRIT (Catalonia), project 2000SGR00025. X.Illa acknowledges financial support from DGICyT.