text:
# Strangeness in Neutron Stars

## 1 INTRODUCTION

The true nature and the internal constitutions of the ultra-dense compact stars known as neutron stars is one of the most fascinating enigma in modern astrophysics . Different models for the equation of state (EOS) of dense matter predict the neutron star maximum mass $`M_{max}`$ to be in the range 1.4 – 2.2 $`M_{}`$ ($`M_{}2\times 10^{33}`$g is the mass of the Sun), and a corresponding central density in the range of 4 – 8 times the normal saturation density ($`\rho _02.8\times 10^{14}`$g/cm3) of nuclear matter. Thus neutron stars are the most likely sites in the universe in which strangeness-bearing matter with a strangeness to baryon ratio $`f_S=S/B1`$ may exist. Strangeness in neutron stars is expected both in a confined form (hyperons and or kaons), or in a deconfined form (strange quark matter). Accordingly different types of “neutron stars” are expected theoretically, as schematically summariezed in Fig. 1.

## 2 NEUTRON STARS OR HYPERON STARS?

In a conservative and oversimplified picture the core of a neutron star is modeled as a uniform fluid of neutron rich nuclear matter in equilibrium with respect to the weak interactions ($`\beta `$–stable nuclear matter). The presence of hyperons in neutron stars was first proposed in 1960 by Ambartsumyan and Saakyan , and since then it has been investigated by many authors. The reason why hyperons are expected in the high dense core of a neutron star is very simple, and it is mainly a consequence of the fermionic nature of nucleons, which makes the nucleon chemical potentials a very rapidly increasing function of density. As soon as the chemical potential of neutrons becomes sufficiently large (see Fig. 2), the most energetic neutrons (i.e. those on the Fermi surface) can decay via the weak interactions into $`\mathrm{\Lambda }`$ hyperons and form a Fermi sea of this new hadronic species with $`\mu _\mathrm{\Lambda }=\mu _n`$. The $`\mathrm{\Sigma }^{}`$ can be produced via the weak process $`e^{}+n\mathrm{\Sigma }^{}+\nu _e`$ when the $`\mathrm{\Sigma }^{}`$ chemical potential fulfill the condition111except from the very initial stage soon after neutron star birth, neutrinos freely escape the star and thus the neutrino chemical potentials have not to be considered in the chemical equilibrium equations. $`\mu _\mathrm{\Sigma }^{}=\mu _n+\mu _e`$. As we can see from the results depicted in Fig. 2, hyperons appear at a relatively moderate density of about 2 times the normal saturation density ($`n_0=0.16`$ fm-3) of nuclear matter. Notice that the $`\mathrm{\Sigma }^{}`$ hyperon appears at a lower density than the $`\mathrm{\Lambda }`$, even though the $`\mathrm{\Sigma }^{}`$ is more massive than the $`\mathrm{\Lambda }`$. This is due to the contribution of the electron chemical potential $`\mu _e`$ to the threshold condition for the $`\mathrm{\Sigma }^{}`$ (i.e. $`M_\mathrm{\Sigma }^{}=\mu _n+\mu _e`$, for free hyperons) and to the fact that $`\mu _e`$ in dense matter is large and can compensate for the mass difference $`M_\mathrm{\Sigma }^{}M_\mathrm{\Lambda }=81.76`$ MeV.

In Fig. 3, we show the profile of a such an hyperon star . As we see the hyperonic matter inner core of the star extend for about 8 km. This radius has to be compared with the total stellar radius $`R11`$ km, and with the thickness of the nuclear matter layer (outer core) which is about 2 km. Thus neutron stars are “giant hypernuclei”  under the influence of gravity and strong interactions.

The influence of hyperons on neutron stars properties has been investigated using different approaches to determine the EOS of hyperonic matter. One of the most popular approaches, to solve this problem, is the relativistic mean field model . Some of the parametrizations of the lagrangian of the theory have tried to reconcile measurd values of neutron star masses with the binding energy of the $`\mathrm{\Lambda }`$ particle in hypernuclei . A different approach is based on the use of local effective potentials to describe the in-medium baryon-baryon (BB) interaction . This method mimic and generalize to the case of hyperonic matter the one based on the Skyrme nuclear interaction in the case of nuclear matter. Here we will report on some recent results based on a third approach, which starts from the basic BB interaction and solve the many-body problen to get the EOS for hyperonic matter. This method is based on an extension of the Brueckener-Bethe-Goldstone (BBG) theory to include hyperonic degrees of freedom . In particular, the study of ref.  focus on the properties of a newborn neutron star, and explore the consequences of neutrino trapping in dense matter on the structural properties and on the early evolution of neutron stars .

In Fig. 4, we show the EOS for $`\beta `$-stable dense matter (i.e. in equilibrium with respect to the weak interactions) obtained by Vidaña et al.  using the Nijmegen soft-core potential (NSC97e) of ref.  to describe the hyperon-nucleon (YN) and hyperon-hyperon (YY) interaction within the Brueckner-Hartree-Fock (BHF) approximation of the extended BBG theory. In ref. , the pure nucleonic contribution to the EOS has been included using a parametrization of the Akmal-Pandharipande EOS , where a semi-phenomenological three-nucleon (NNN) interaction of the Urbana type is added to the nuclear hamiltonian to reproduce the empirical saturation point of nuclear matter.

As expected, the presence of hyperons makes the EOS much softer with respect to the pure nucleonic case. The softening of the EOS caused by the presence of hyperons has important consequences on many macroscopic properties of the star: the maximum stellar mass is reduced by $`\mathrm{\Delta }M_{max}`$ 0.5 – 0.8 $`M_{}`$, and the corresponding central density is increased. Also, hyperon stars are more compact (i.e. they have a smaller radius) with respect to pure nucleonic neutron stars. This is illustrated in Fig. 5, where we show the mass radius-relation for traditional neutron stars and for hyperon stars obtained with the microscopic EOS of ref.  (left panel) and with the relativistic mean field EOS (GM3 model) given in ref. . The results depicted in Fig. 5 clearly demonstrate that to neglect hyperons leads to an overstimate of the maximum mass of neutron stars.

It is important to notice the “low” value of the stellar maximum mass, predicted within the approach of ref. , which is in contrast with some precise determination of neutron star masses . For example in the case of the neutron star associated to the pulsar PSR1913+16, the measured stellar mass is

$$M_{PSR1913+16}=(1.4411\pm 0.0007)M_{}.$$
(1)
The prediction of a value for $`M_{max}`$ below the measured neutron star masses is a common feature of all the present microscopic EOS of hyperonic matter based on G-matrix BHF calculations . For example, the authors of ref. , in case of the Argonne $`v_{18}`$ NN interaction, found $`M_{max}=2.00M_{}`$, a corresponding radius of $`R=10.54`$ km and a central density $`\rho _c=1.11`$ fm-3 for neutron stars with a pure nucleonic core. When hyperons are considered as possible stellar constituents, they found  $`M_{max}=1.22M_{}`$, a corresponding radius of $`R=10.46`$ km and a central density $`\rho _c=1.25`$ fm-3. Therefore the current equations of state for hyperonic matter, deduced from microscopic G-matrix BHF calculations, are “too soft” to explain observed neutron star masses.

Clearly, one should try to trace the origin of this problem back to the underlying YN and YY two body interactions or to the possible repulsive three-body baryonic forces involving one or more hyperons, not included in the work of ref. . Presently this is a subject of very active research by people working in this field. Therefore, the use of microscopic EOSs of hyperonic matter in the contest of neutron star physics is of fundamental importance for our understanding of the strong interactions involving hyperons, and to learn how these interactions behave in dense many-body systems.

## 3 KAON CONDENSATION IN NEUTRON STARS

The inner core of neutron stars could also contain a Bose-Einstein condensate of negative kaons . As the density of stellar matter is increased, the $`K^{}`$ energy is lowered by the attractive vector mean field originating from dense nucleonic matter. When the $`K^{}`$ energy becomes smaller than the electron chemical potential $`\mu _e`$ (which is an increasing function of density) the strangeness changing process $`e^{}K^{}+\nu `$ becomes possible. The critical density for this process has been calculated to be in the range 2.5 – 5.0 $`n_0`$ .

Due to the lack of space, we do no have the possibilty to discuss the many relevant implications that kaon condensation has for the structure and the evolution of neutron stars. We refer the reader to the original literature on the subject (see e.g.  and references therein quoted).

## 4 HYBRID STARS

The core of the more massive neutron stars is one of the best candidates in the Universe where a phase transition from hadronic matter to a deconfined quark phase should occur. The quark-deconfinement phase transition proceeds through a mixed phase over a finite range of pressures and densities . At the onset of the mixed phase, quark matter droplets form a Coulomb lattice embedded in a sea of hadrons and in a roughly uniform sea of electrons and muons. As the pressure increases various geometrical shapes (rods, plates) of the less abundant phase immersed in the dominant one are expected. Finally the system turns into uniform quark matter at the highest pressure of the mixed phase . Compact stars which possess a quark matter core, either as a mixed phase of deconfined quarks and hadrons, or as a pure quark matter phase, are called Hybrid Stars .

Many possible astrophysical signals for the appearence of a quark core in neutron stars have been proposed in the last few years (see  and references therein quoted). Particularly, pulse timing properties of pulsars have attracted much attention since they are a manifestation of the rotational properties of the associated neutron star. The onset of quark-deconfinement in the core of the star, will cause a change in the stellar moment of inertia . This change will produce a peculiar evolution of the stellar rotational period ($`P=2\pi /\mathrm{\Omega }`$) which will cause large deviations of the so called pulsar braking index $`n(\mathrm{\Omega })=(\mathrm{\Omega }\ddot{\mathrm{\Omega }}/\dot{\mathrm{\Omega }}^2)`$ from the canonical value $`n=3`$, derived within the magnetic dipole model for pulsars and assuming a constant moment of inertia for the star. The possible measurement of a value of the braking index very different from the canonical value (i.e. $`|n|>>3`$) has been proposed  as a signature for the occurrence of the quark-deconfinement phase transition in a neutron star. However, it must be stressed that a large value of the braking index could also results from the pulsar magnetic field decay and/or alignment of the magnetic axis with the rotation axis .

## 5 STRANGE STARS

The possible existence of a new class of compact stars completely made of deconfined u,d,s quark matter (strange quark matter (SQM)) is one of the consequences of an hypothesis  formulated by A.R. Bodmer in 1971 and revived by E. Witten in 1984. These stars are usually called Strange Stars. According to the Bodmer-Witten hypothesis SQM could be the true ground state of matter. In other words, at zero temperature and pressure, the energy per baryon of SQM could be less than the energy per baryon of 56Fe, which is the most tightly bound nucleus in nature. The strange matter hypothesis does not conflict with the existence of atomic nuclei as conglomerates of nucleons, or with the stability of “ordinary” matter . Thus strange stars may exist in the universe.

One of the most likely strange star candidate is the compact object in the transient X-ray burst source SAX J1808.4-3658 (ref. ). This X-ray source was discovered in 1996 by the BeppoSAX satellite. Two bright type-I X-ray bursts were detected, each lasting less than 30 seconds. Analysis of the bursts in SAX J1808.4-3658 indicates that it has a peak X-ray luminosity of $`6\times 10^{36}`$erg/s in its bright state, and a X-ray luminosity lower than $`10^{35}`$erg/s in quiescence. SAX J1808.4-3658 is a X-ray millisecond pulsar with a pulsation period of 2.49 ms, also it is a member of a binary stellar system with orbital period of two hours. Using the observational data collected by the Rossi X-ray Timing Explorer during the the 1998 April-May outburst, Li et al.  have obtained an upper limit for the compact star radiud as a function of the unknown stellar mass. Comparing this observational mass-radius (M-R) relation of SAX J1808.4-3658 with the theoretical M-R realtions for traditional neutron stars, hyperon stars, stars with kaon condensation, and strange stars Li et al.  (see their Fig. 1) argue that a strange star model is more consistent with SAX J1808.4-3658, and suggest that it could be a strange star.

SAX J1808.4-3658 is not the only LMXBs which could harbour a strange star. Recent studies have shown that the compact stars associated with the X-ray burster 4U 1820-30 (ref.), the bursting X-ray pulsar GRO J1744-28 (ref.), the X-ray pulsar Her X-1 (ref.), the kHz QPOs source 4U 1728-34 (ref.), are likely strange star candidates. Recently, it has been suggested that the isolated compact star RX J1856.5-3754 (ref. ) could be a strange star.

## 6 QUARK-DECONFINEMENT PHASE TRANSITION IN NEUTRON
STARS AND GAMMA-RAY BURSTS

Gamma Ray Bursts (GRBs) are one of the most violent and mysterious phenomena in the universe (see e.g. ref. for a general introduction on this subject). During the last ten years two satellites, the Compton Gamma Ray Observatory (CGRO) and BeppoSAX, have revolutionized our understanding of GRBs. The Burst And Transient Source Experiment (BATSE) on board of the CGRO has demonstrated that GRBs originate at cosmological distances. The BeppoSAX discovered the X-ray afterglow. This has permitted to determine the position of some GRBs, to identify the host galaxy and, in a number of cases, to measure the red-shift. If the energy is emitted isotropically, the measured fluence of the bursts implies an energy of the order of $`10^{53}`$ erg. However, there is now compelling evidence that the $`\gamma `$-ray emission is not isotropic, but displays a jet-like geometry. In this case the GRB energy is of the order of $`10^{51}`$ erg .

Many cosmological models for the energy source of GRBs have been proposed. Presently one of the most popular is the so-called ”collapsar”, or ”hypernova” model. Alternative models are the merging of two neutron stars (or a neutron star and a black hole) in a binary system, or the accretion of matter into a black hole. The present report is not the appropriate place to discuss the various merits and drawbacks of the many theoretical models for GRBs. In the following, we will mention some recent research which try to make a connection between GRBs and quark-deconfinement phase transition.

A possible central engine for GRBs is the conversion of a pure hadronic compact star to a strange star. The stellar conversion is triggered by the formation of a SQM drop in the center of the hadronic star. This idea was proposed long time ago by Alcook et al. . Recently detailed calculations, based on different realistic models for the equation of state of neutron star matter and SQM, have been performed by the authors of ref.. They showed that the total amount of energy liberated in the conversion is in the range (1 – 4)$`\times 10^{53}`$erg. This energy will be mainly taken away by the neutrinos produced during the quark-deconfinement phase transition. If the efficiency of the conversion of neutrinos to $`\gamma `$ is of the order of a few percent , then the birth of a strange star from a neutron star could be the energy source for GRBs.

A mounting number of observational data suggest a clear connection between supernova (SN) explosions and GRBs . Particularly, in the case of the gamma ray burst of July 5, 1999 (GRB990705) and in the case of GRB011211, it has been possible to estimate the time delay between the two events. For GRB990705 the supernova explosion is evaluated to have occurred a few years before the GRB , while for GRB011211 about four days before the burst .

The scenario which emerges from these findings is the following two-stage scenario: (i) the first event is the supernova explosion which forms a compact stellar remnant, i.e. a neutron star (NS); (ii) the second catastrophic event is associated with the NS and it is the energy source for the observed GRB. These new observational data, and the scenario outlined above, poses severe problems for most of the current theoretical models for the central energy source of GRBs. The main difficulty of all these models is to give an answer to the following questions: what is the origin of the second “explosion”? How to explain the long time delay between the two events?

In the so-called supranova model  for GRBs the second catastrophic event is the collapse to a black hole of a supramassive neutron star, i.e. a fast rotating NS with a baryonic mass $`M_B`$ above the maximum baryonic mass $`M_{B,max}`$ for non-rotating configurations. In this model, the time delay between the SN explosion and the GRB is equal to the time needed by the fast rotating newly formed neutron star to get rid of angular momentum and to reach the limit for instability against quasi-radial modes where the collapse to a black hole occurs . The supranova model needs a fine tuning in the initial spin period $`P_{in}`$ and baryonic stellar mass $`M_{B,in}`$ to produce a supramassive neutron star that can be stabilized by rotation up to a few years. For example, if $`P_{in}1.5`$ ms, then the newborn supramassive neutron star must be formed within $`0.03M_{}`$ above $`M_{B,max}`$ .

In a very recent paper, Berezhiani et al.  (see also ref. ) have proposed a new model to explain the SN–GRB association and in particular the long time delay inferred for GRB990705 and GRB011211. In the model of ref., the second explosion is related to the conversion from a metastable purely Hadronic Star (neutron star or hyperon star) into a more compact star in which deconfined quark matter is present (i.e. a hybrid star or a strange star). The new and crucial idea in the work of ref.  with respect to previous work , is the metastability of the hadronic star due to the existence of a non-vanishing surface tension at the interface separating hadronic matter from quark matter. The mean-life time of the metastable hadronic star can then be connected to the delay between the SN explosion and the GRB. The nucleation time (i.e. the time to form a critical-size drop of quark matter) can be extremely long if the mass of the star is small. Via mass accretion the nucleation time can be dramatically reduced and the star is finally converted from the metastable into the stable configuration . A huge amount of energy, of the order of 1052–1053 erg, is released during the conversion process and can produce a powerful gamma ray burst. Within the model proposed by Berezhiani et al.  is is possible to have different time delays between the two events since the mean-life time of the metastable hadronic star depends on the value of the stellar central pressure. Thus the model of ref.  is able to interpret a time delay of a few years (as observed in GRB990705 ), of a few days (as in the case of GRB011211 ), or the nearly simultaneity of the two events (as in the case of SN2003dh and GRB030329 ).

## Acknowledgments

I am greatly indebted to I. Vidaña for providing me some unpublished results from his Ph.D. thesis and for the great help in preparing most of the figures of the present work.


format: text

text:
# Using virtual processors for SPMD parallel programs

## 1 Introduction

The processors virtualization technology permits to split a real physical processor into two virtual chips, so that the operating system, as MS Windows or Linux, of a computer can use the virtual processors as two real chips. Example of such technology is Intel’s Hyper Threading . The hardware can so be considered as a symmetric multi-processor machine and the software can use it as a true parallel environment.

In this work I show some results obtained with parallel computations using Matlab  programs on Intel technology. The physical and logical characteristics of the used machine are presented in the following tables:

| Hardware |  |
| --- | --- |
| Type | HP Compaq ProLiant DL380 |
| Processors | 2 Intel Xeon 2.40 GHz |
| Ram | 1.5 GB |
| Storage | 6 SCSI disks 36.5 GB - Raid 5 |

| Software |  |
| --- | --- |
| Operating System | MS Windows Server 2000 |
| Matlab | v. 5.3 |

The Matlab programs used for these experiments was based on cycles of floating-point computations.

## 2 The parallel Matlab environment

The package Matlab has not a native support for parallel elaboration and multithreading . Yet, there are some extensions, as tools and libraries , for the use of a parallel environment on multi-processors hardware. With the biprocessor machine I have preferred the method of splitting a given computation on multiple instances of the runtime Matlab program. A single master instance starts the slave copies and assigns to each of them the same set of instructions on different sets of data. Hence I have simulated a SPMD computation on the machine.

In this way the parallel environment is simple, because there is not need of external libraries or calls to interfaces, and flexible, because to a single slave copy it can be assigned a set of different instructions for realizing a MPMD computation.

With this method the exchange of messages among independent processes is a problem. The only way to communicate from one Matlab copy to another is the use of shared files. In a second type of experiments I show that this method is not critical for the time execution if one uses fast mass-storage as SCSI or FiberChannel systems.

### 2.1 The SPMD programs

In the experiments I have defined a master Matlab function which writes to a shared file system the .m scripts to be executed by slaves Matlab copies. These copies are launched in background mode for the parallel execution. The master program controls the end of the computations using a simple set of lock-files. The slaves finish their work, save on files the results and cancel the own lock-file. The master reads the sets of data from these files for other possible computations. Now I describe the principal code of the program.

This is the declaration of the function masterf. The lockarray variable is an array for testing the presence of the lock-files during the slaves computation. The finalres is an array for the collection of the partial results from slaves. The string computing is the mathematical expression to use in the computation.

function \[elapsedtime,totaltime,executiontime\]=masterf(nproc,maxvalue,step,computing)
%
% MASTERF: master function for parallel background computation.
%
% sintax:
% \[elapsedtime,totaltime,executiontime\]=masterf(nproc,maxvalue,step,computing)
%
% input parameters:
%
% nproc = number of processes;
% maxvalue = sup-limitation of the data-array to process; the inf-limitation is 0;
% step = difference from two consecutive numbers in the data-array;
% computing = the string of the mathematical expression to compute;
%
% output parameters:
%
% elapsedtime = total elapsed time to complete the execution of the computation;
% totaltime = sum of the single slaves CPU-time to complete the single computation;
% executiontime = single slaves CPU-time to complete the assigned computation;
ostype=computer;
tottime=0.;
lockarray=0:nproc-1;
numbervalues=maxvalue;
computingstring=\[’ ’ computing\];
finalres=\[\];

After the assignment of the own value to variable workdir, working directory of Matlab, a cycle writes on storage the slaves lock-files.

for i=0:nproc-1
filelock = strcat(workdir,’filelock’,int2str(i));
fid=fopen(filelock,’wr’);
fwrite(fid,”);
fclose(fid);
end

In the next fragment of program, the master sets the commands for the writing of an appropriate Matlab .m script for every slave process. Such script contains the instruction for determining the CPU-time spent on calculus, the expression of the mathematical computation, the instruction to save on storage the data computed and the CPU-time, finally the instruction to delete the lock-file.

for i=0:nproc-1
if (i==0) middlestep=0; else middlestep=1; end
infdata=i\*(numbervalues/nproc) + middlestep\*step;
supdata=(i+1)\*(numbervalues/nproc);
fileworker = strcat(workdir,’fileworker’,int2str(i),’.m’);
commandworkertmp = …
strcat(’x=’,num2str(infdata),’:’,num2str(step),’:’,num2str(supdata),…
’; t1=cputime; ’,computingstring,…
’; t2=cputime-t1; save out’,int2str(i));
commandworker = \[’cd ’ workdir ’; ’ commandworkertmp …
’ y t2; ’ ’delete filelock’int2str(i) ’; exit;’\];
fid = fopen(fileworker,’wt’);
fwrite(fid,commandworker);
fclose(fid);
end

After the instructions for determining the CPU-time and the elapsed-time (tic) spent by the master program, a cycle launches the slaves Matlab runtimes. In the case of Windows operating system, the startcommand string is ”start”, an OS command for the background running of an executable program, and the osstring string is ”dos”. In the case of Unix-like operating system, the string are ”sh” and ”unix” respectively. Each slave executes immediately the fileworker script, as shown by the Matlab ”-r” parameter.

t1 = cputime;
tic;
for i=0:nproc-1
fileworker = strcat(’fileworker’,int2str(i));
commandrun = \[startcommand ’ matlab -minimize -r ’ fileworker\];
eval(strcat(\[osstring,’(’,””,commandrun,””,’);’\]));
end

In the next fragment of code the master program executes a cycle for determining the end of slaves computations. It controls if the lockarray variable has some process’s rank non negative. In this case, it attempts to open the relative lock-file; if the file still exists, the master closes it, else the lockarray process position is set to -1. The pause instruction can be useful for avoiding an excessive frequency, hence an high cpu-time consuming, in the ”while” cycle.

lockarraytmp=find(lockarray $`>`$ -1);
while (length(lockarraytmp) $`>`$ 0)
pause(.1);
for i=lockarraytmp
fid = fopen(strcat(’filelock’,int2str(i-1)),’r’);
if (fid $`<`$ 0)
lockarray(i) = -1;
else
fclose(fid);
end
end
lockarraytmp=find(lockarray $`>`$ -1);
end

At the end, the master reads the partial slaves computation outputs and stores them in an array. At this point the master cpu-time and elapsed time are registered too. The total execution time is defined as sum of the slaves computation cpu-time, and is useful for comparison with the execution time in the case $`\text{nproc}=1`$. The single slave execution time is defined as the arithmetic mean of all the partial execution times.

for i=0:nproc-1
partialres = load(strcat(’out’,int2str(i)));
finalres = \[finalres partialres\];
end
elapsedtime = toc;
totaltime = cputime - t1;
for i=0:nproc-1
fps = load(strcat(’out’,int2str(i)));
tottime = tottime + fps.t2;
executiontime = tottime/nproc;
end

## 3 Tests and results

For the tests I have used the following values for the masterf parameters:
nproc: from 1 to 8;
maxvalue: m * 10000, where m = 1, 2, 3;
step: 0.001;
computing: $`y=5432.060708\mathrm{cos}((\mathrm{sin}(x^{9.876}))^{1.2345})`$.

I have also tested the program without the slaves saving of partial computations results and their final master load, for determining the influence of the I/O storage operations on the times of execution.

The following are the results obtained with four tests for every type of experiment. The values are arithmetic means approximated to two decimals and they are expressed in seconds; the numbers from 1 to 8 are the value of nproc, while ”m” is the factor parameter in the maxvalue expression.
I have not reported the elapsed-times, because they weren’t different from the cpu-times registered, probably due to the fact that, during the experiments, the server was dedicated only to the computations.

Tables of results.
All the values are expressed in seconds.

1.a Medium execution cpu-times for process, no data storage:

| m | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | 39.32 | 20.89 | 14.60 | 11.49 | 11.86 | 8.73 | 9.15 | 7.29 |
| 2 | 77.56 | 40.83 | 29.20 | 23.41 | 22.31 | 23.19 | 17.93 | 19.02 |
| 3 | 137.75 | 69.70 | 51.67 | 35.28 | 34.91 | 36.98 | 32.46 | 34.86 |

2.a Total execution cpu-times, no data storage:

| m | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | 41.01 | 22.11 | 16.19 | 13.67 | 18.42 | 17.10 | 18.17 | 16.20 |
| 2 | 78.40 | 41.89 | 30.74 | 24.81 | 26.78 | 34.69 | 33.68 | 32.95 |
| 3 | 139.05 | 75.66 | 59.38 | 38.83 | 40.28 | 49.69 | 48.98 | 48.27 |

1.b Medium execution cpu-times for process, with data storage:

| m | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | 42.78 | 20.59 | 15.99 | 14.61 | 11.05 | 8.89 | 11.56 | 10.94 |
| 2 | 99.93 | 52.66 | 38.59 | 25.22 | 20.71 | 18.60 | 18.55 | 18.70 |
| 3 | 151.03 | 80.49 | 57.33 | 39.65 | 28.83 | 36.03 | 44.16 | 44.76 |

2.b Total execution cpu-times, with data storage:

| m | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1 | 49.55 | 28.70 | 26.69 | 27.26 | 17.98 | 18.92 | 17.78 | 17.76 |
| 2 | 131.66 | 68.50 | 54.33 | 40.53 | 38.00 | 38.37 | 36.76 | 40.69 |
| 3 | 201.65 | 102.90 | 75.03 | 58.94 | 60.97 | 66.52 | 67.92 | 67.83 |

## 4 Analysis of results

From the results of the previous section, I deduce the following observations:

1. In table 1.a the gain in execution speed is good from 1 to 4 processes, while from 5 to 8, and in particular in the case m=3, the gain is low; this fact can be due to the excess load on the dual-processor machine when nproc $`>`$ 4;
2. In the same table the speedup  is quasi-linear from 1 to 4, hence the algorithm and the parallel environment used are an optimized SPMD implementation if one is interested only on pure computation time;
3. In table 2.a the value nproc=4 gives the best performance as total execution time; hence, if one is interested on the time spent by the master to verify when the slaves finish their job, the four virtual processors guaranteed by Hyper Threading technology show the best efficiency in the case of four running processes;
4. With exclusion of the case nproc=1, when the master must verify only a single process, in the case m=3 the difference between the times rows of 2.a and 1.a is smallest in the case nproc=4; hence in this case too the total time registered by the master is optimized respect to the slaves execution time;
5. In the case of data storage, the preceding conclusions are not so clear, probably due to the fact that the reading-writing of small files is in general optimized by the RAID 5 technology on a multi-disks system; this fact seems to be confirmed by the great difference, 50 seconds, from the times registered in the case m=3 and nproc=1 in tables .b;
6. In table 2.b one can note that for m=3 the best result is in the case nproc=4;
7. In tables .b, in the case m=3 the difference between the times rows is smallest in the case nproc=4 too;
8. From tables 2.b, one can note that for m=1 the best speedup (8 processes) respect to the case nproc=1 is 2.78, while for m=3 the best (4 processes) is 3.42; hence the virtual processors seem to have a better performance with a large amount of data to compute, probably due to the fact that in this case the parallel computations have a greater relevance respect to the physical operations on storage system.

### 4.1 Conclusions

From the previous facts one can deduce that a virtual processors technology as Hyper Threading can be a good choice for running SPMD programs in the case that

* the number of parallel processes is equal to the number of virtual processors;
* the data to be computed have a large amount, particularly when their distribution among processes and the merging of final results are based on files stored on fast storage system.

## 5 Acknowledgements

I wish to thank Ernesto Montagnoli, Mario Bubola and Filippo Rossetti, chief and members of Infrastructures and Systems Area at ICT Department, Riello Group, for the permission to execute this work on the ProLiant server and for technical assistence.


text:
# Start-to-End Simulations of SASE FEL at the TESLA Test Facility, Phase 1

## 1 Introduction

The vacuum ultra-violet (VUV) free-electron laser (FEL) at the TESLA Test Facility (TTF), Phase 1 demonstrated saturation in the wavelength range 80-120 nm based on the self-amplified spontaneous emission (SASE) principle. Analysis of experimental data for the radiation properties have led us to the unique conclusion that SASE FEL produced ultra-short radiation pulse (FWHM duration 30-100 fs) with GW level of peak power. However the measured properties of the radiation were in strong disagreement with project parameters for the electron bunch : we expected a peak beam current of 500 A, an rms bunch duration of 1 ps, and an energy spread of 0.5 MeV. Such electron beam parameters would result in operation of the FEL in the saturation regime only for normalized emittance less or about 2 mm-mrad, while emittance measurements gave values of approximately 4 mm-mrad for a bunch charge of 1 nC . Also, the radiation pulse was by one order of magnitude shorter than the value expected for the aforementioned project parameters.

Facing this evident disagreement, we did an attempt to reconstruct parameters of the lasing part of the electron bunch using measured properties of the radiation. Actually, radiation measurements were very reliable and accurate, and in combination with the FEL theory we can infer a lot about the properties of the part of electron bunch that produces the radiation. In the FEL simulations lasing part of the bunch was approximated by a Gaussian. A set of parameters for this lasing part of the electron bunch leading to simulation results, consistent with the radiation measurements, were: a FWHM bunch duration of 120 fs, a peak current of 1.3 kA, an rms energy spread of 100 keV, and a rms normalized transverse emittance of 6 mm-mrad. These values, as used in the FEL-simulation, gave reasonable agreement for the average energy in the radiation pulse and statistical distributions of the fluctuations of the radiation energy. However, there were two visible disagreements with experimentally measured radiation characteristics. First, there was a noticeable difference in the shape of angular distribution of the radiation intensity. Second, the measured averaged spectrum width was visibly wider with respect to the simulated one and the spikes in the single-shot spectra were larger in the experiment than the one simulated.

It is worth mentioning that soon after obtaining the saturation, attempts for more detailed electron beam measurements were undertaken. The first one was measurement of the slice energy spread before compression  which gave the value of about 5 keV. This measurement clearly indicated that the value of peak current after compression should be well above 1.3 kA. The second experiment aimed at direct time-domain measurement with a streak camera  of the bunch shape and peak current. The results of these measurements gave $`650\pm 100`$ fs for rms pulse duration, and 0.6 kA peak current at 3 nC. These results were consistent with another techniques based on a tomographic reconstruction of the longitudinal phase space. Both of these measurement were in fact resolution limited: for instance the temporal resolution of the streak camera in the experiment reported in  was 200 fs. It was therefore impossible to resolve any fine structure on the bunch charge density that have characteristic width below $``$200 fs. Hence it was impossible to resolve the lasing fraction of the bunch since it was below 200 fs as we mentioned above. We should also note that there was no possibility to measure slice emittance and energy spread at TTF.

Thus there was no complete quantitative description of the TTF SASE FEL operation. To get a better understanding of the TTF SASE FEL operation we undertook the full physics simulation of the TTF FEL, Phase 1. This study aimed to trace the evolution of the electron beam from the photo-cathode to the undulator entrance. Then use the thereby produced electron distribution to calculate the radiation produced by this electron bunch while passing in the undulator. In our simulations we tried (to our best knowledge) to reproduce the main parameters of the machine from the FEL run in September 2001. In some cases, due to the lack of reliable informations, we had to make simplified assumptions.

Presently there is no universal particle tracking code capable of the evolution of the electron beam through the accelerator and bunch compression chain. For the simulation reported hereafter, the program ASTRA , which takes into account space charge, but does not calculate the beam dynamics in the bends, was used in straight transport sections (gun, capture cavity, accelerating modules, drift spaces). The beam dynamics in the bunch compressor was independently simulated with the codes elegant  and CSRtrack  taking into account the effects related to coherent synchrotron radiation (CSR). Finally, the FEL process was simulated by three-dimensional, time-dependent code FAST .

## 2 Facility description

The description of the TTF accelerator (see Fig. 1), operating under standard lasing conditions, can be found in and references therein. The RF gun consists of an L-band cavity (1+1/2 cell) incorporating a CsTe2 photo-cathode illuminated by a UV laser with a (Gaussian) pulse distribution of 7-8 ps rms. The electron bunch with the charge 2.8 nC and energy about 4 MeV is extracted from the gun (nominal laser launch phase is 40 deg) and is then accelerated in the booster cavity up to 16 MeV. The phase of the booster cavity is normally chosen such that the total (correlated) energy spread is minimized. Passing a rather long drift, the beam is then injected into the superconducting TESLA module (ACC1) where it is accelerated up to 135 MeV off-crest to impart the proper correlated energy spread for subsequent compression in the following four-bend magnetic chicane (BC2). Without compression the bunch is rather long (about 3.5 mm rms - longer than the laser pulse due to Coulomb repulsing in the injector). Thus, at the nominal compression phase (10 deg off-crest) the bunch accumulates RF curvature leading to a ”banana” shape on the longitudinal phase plane after compression. The resulting bunch shape in time domain constitutes a short high-current leading peak and a low-current long tail. After compression the beam passes a 5 m long drift, the second TESLA module ACC2 (being accelerated up to 248 MeV), and another drift (about 20 m) that includes a collimation section and a transverse matching section. Finally the beam enters the undulator consisting of three 4.5 m long modules where a short SASE pulse (wavelength below 100 nm) is produced by the leading peak of the bunch. The electron beam is separated from the photon beam thanks to the spectrometer dipole, and bent in a beam dump while the photon beam goes to the photon diagnostic area.

## 3 Simulations of the beam dynamics in the accelerator

The initial part of the accelerator, from the photo-cathode to BC2 entrance, was simulated with Astra , a program that includes space charge field using a cylindrical symmetric grid algorithm. The beam was then tracked through the bunch compressor with elegant  that includes a simplified model (based on a line charge approximation) of the CSR wake. Downstream of BC2, Astra was again used up to the undulator entrance because the space charge induced-effects are significant for the strongly compressed part of the bunch. SASE FEL process in the undulator was simulated with three-dimensional time-dependent code FAST . In order to check that we did not miss any important CSR-related effects in the bunch compressor, we performed an independent simulations of the bunch compressor with a newly developed code CSRtrack . This latter code incorporates a two-dimensional self-consistent model of the beam dynamics.

The Astra simulation from the cathode to the bunch compressor was performed using $`2\times 10^5`$ macro-particles. The main parameters of the electron beam, calculated for various longitudinal slices along the bunch, at the entrance to ACC1 are presented in Fig. 2. The chosen bin size is such that the energy chirp does not contribute to the slice energy spread. The latter parameter is very important since it defines the width of the leading peak and the peak current after compression. One can see that it varies along the bunch, taking the value of about 3.5 keV at the location $`s=0`$ (this part of the bunch is then put into a local full compression in BC2). The resulting values for the local energy spread in Fig. 2 are in agreement with the measurements . The main source of the build-up of the local energy spread (as we see it from this simulation) is a transverse variation of the longitudinal space charge field in the injector. Since it is a coherent effect (a particle’s energy deviation is correlated with its position in the bunch, in particular, with its transverse offset within a given slice), the frequently used notions of ”uncorrelated” or ”incoherent” energy spread are not adequate here. As we will see, the correlations can be important when we compress the beam.

At the entrance of BC2 the output distribution of macro-particles is converted into the input distribution for elegant. It is then tracked through the bunch compressor ($`R_{56}23`$ cm) without CSR wake included. The resulting distribution on longitudinal phase plane and slice parameters for the leading peak are presented in Fig. 4. One can see that the peak current and the width of the current spike are in a reasonable agreement with analytical estimate, presented in Appendix. The difference can be explained by the non-Gaussian energy distribution in a slice before compression and by the fact that the beam entered ACC1 with the energy 16 MeV and some energy chirp, etc. Also, important correlations in particle distribution are not considered in Appendix. What requires some explanations is a behavior of the slice emittance in Fig. 4. As we have already mentioned, the main effect, responsible for a slice energy spread in the injector, is the transverse variation of the longitudinal space charge field. An energy offset of a particle in a slice (with respect to on-axis particle) is correlated with its transverse offset (and, since particles are moving transversely, with a betatron amplitude). A sign of the energy offset depends on the position of the slice in the bunch. Indeed, the head of the bunch is accelerated due to space charge field, so that for a given slice a particle with a finite betatron amplitude gets less energy than an on-axis particle. The situation is opposite for the slices which are on the left from a minimum energy spread in Fig. 2: the larger transverse offset, the larger positive energy deviation. During compression particles with higher energies (in our case, with larger betatron amplitudes) move forward, this explains the large slice emittance values for the slices in the leading front of the spike. On the other hand, due to this effect a slice with the maximal current has less ”bad” particles, so that the emittance there is visibly less than that at $`\mathrm{s}0`$ in Fig. 2. We also tried to put the head of the bunch (say $`\mathrm{s}5`$ mm in Fig. 2) into a local full compression. The result was that after compression emittance was smoothly decreasing in the leading front of the spike. It is worth noting that when we prepared a 6-D Gaussian bunch with no correlations (at the entrance to ACC1) and track it through BC2, the slice emittance after compression was the same in all slices and was equal to its value before compression.

Our next step was to track the same particle distribution through the bunch compressor, taking into account CSR wake in a simplified model used in elegant. The beam parameters behind BC2 are presented in Fig. 4. Note that due to CSR-induced effects the peak current slightly decreased, by less than 10%. Slice emittance (in horizontal plane) in the slice with maximal current increased by some 50%, and the local energy spread by almost a factor of 2.

It looks surprising that the peak current is almost unchanged in the presence of the CSR wake. Indeed without including CSR, the final shape of the bunch forms in the end of the third - begin of the fourth dipole of BC2. If one does a naive estimate for the energy kicks due to CSR (on the way to the end of the fourth dipole) for such a narrow peak with so high current, and applies $`R_{56}`$ to the end of compressor, then one finds out that the distribution of current should be strongly disturbed. The reason why this does not happen can be explained as follows. For our range of parameters, one can not neglect coupling between transverse and longitudinal phase spaces in the bunch compressor (the importance of this effect was pointed out in studies of CSR microbunching instability ), described by linear transfer matrix elements $`R_{51}`$ and $`R_{52}`$ (the net effect through the whole compressor is zero to first order, and higher order terms are negligible). This coupling makes the leading peak effectively much longer and the current much lower than in the case of zero emittance (the spike gets cleaned up only at the very end of the fourth dipole). Thus, the energy kicks due to CSR are strongly suppressed, and the current spike survives.

In order to confirm this result and to check that we did not miss any important effect, using a simplified CSR model in elegant, we performed alternative simulations for the bunch compressor with the ”first principles” code CSRtrack. For technical reasons, the incoming distribution of macro-particles has to be simplified in those simulations (for instance, the above mentioned correlations were neglected), so that slice parameters after compression are not exactly reproduced even without CSR. Nevertheless, the main results of simulations with elegant were confirmed: the peak current decreases by less than 10%, and the slice emittance growth is in the range of several tens of per cent.

¿From BC2 exit to the undulator entrance the tracking was done with Astra. The reason is that a simple estimate predicts very strong longitudinal space charge effect in the leading peak. Indeed, for a Gaussian bunch with an rms length $`\sigma _z`$ and a peak current $`I`$ the change of the peak-to-peak energy chirp $`\mathrm{\Delta }\gamma `$ (in units of the rest energy) in a drift can be estimated as

$$\frac{d(\mathrm{\Delta }\gamma )}{dz}2.4\frac{I}{I_\mathrm{A}}\frac{\mathrm{ln}(\gamma \sigma _z/\sigma _{})}{\sigma _z\gamma ^2}$$

where $`I_\mathrm{A}=17`$kA is the Alfven current, $`\gamma `$ the relativistic factor and $`\sigma _{}`$ the rms transverse size of the beam. This formula holds when $`\sigma _z\gamma \sigma _{}`$. The estimate shows that for the leading peak the energy chirp should be in the range of several MeV.

To reduce numerical calculation effort, we did not track the entire distribution since we were anyway not interested in the parameters of the long low-current tail. So, we cut the tail away and tracked particles in the head of the bunch (typically this part was a few hundred $`\mu `$m long in our simulations). The parameters of the front part of the bunch at the undulator entrance are shown in Fig. 5. One can notice a big energy chirp due to the space charge within the current spike. Note also that due to Coulomb repulsing the spike gets wider, and the peak current decreases by approximately 20%. A change of the local energy spread is due to transverse variations of the longitudinal space charge field. As one would expect, the local minimum of the energy spread is close to the position of maximal current (where derivative of the current is zero).

## 4 Simulation of the SASE FEL process

The simulations SASE FEL process were performed with the three-dimensional time-dependent FEL code FAST . Before describing the results of these simulations we should discuss a difficulty, connected with the data transfer from a beam dynamics simulation code to an FEL simulation code. While the data exchange between beam dynamics simulation codes (for instance Astra $``$ elegant and elegant $``$ Astra) is straightforward and technically simple, a direct loading, for instance, of Astra output distribution as FAST input distribution is impossible. The reasons for this are: completely different time scales of the processes, and a necessity to avoid an artificial noise of macro-particles in FEL spectral range (in addition, in case of a SASE FEL simulation it is also necessary to correctly simulate a real shot noise in the beam - an input signal for SASE FEL).

The standard way of preparing input data for an FEL code is as follows: a macro-particle distribution at the undulator entrance is cut into longitudinal slices. A mean energy, rms energy spread, current, rms emittances etc. are calculated for each slice. Then regular distributions (normally - Gaussian) in each slice are used in FEL code for energy and both transverse phase spaces, having the same mean and rms values as input distribution. Clearly, since only mean and rms values are extracted from the original distributions, an essential information can be lost if the distributions are strongly non-Gaussian. For energy distribution it is sufficient to know only mean and rms values if the rms value is visibly less than FEL parameter $`\rho `$ , which is often the case. For transverse phase space, however, such a simplified procedure may not be satisfactory, especially when the rms emittance is strongly influenced by ”halo” particles. More adequate procedure could be a two-dimensional Gaussian fit (with the least-square method) for each longitudinal slice of each phase plane. Indeed, not only first and second, but also higher-order momenta of the distribution contribute to the fit result. A disadvantage of this procedure is the requirement to have a lot of macro-particles in each slice in order to do a reliable fit. To avoid this difficulty we did the following. We used Gaussian fit in $`x`$ and $`y`$-phase spaces for the part of the bunch shown in Fig. 5, without cutting it into slices, and got new values for rms emittances. Then the slice emittance curves were rescaled in accordance with the fit result. Finally, after some smoothing we got the input parameters for FAST which are shown in Fig. 6. In the simulation we assumed that all slices are perfectly matched to the undulator, and all the centroids are on the ideal orbit.

We performed 300 statistically independent runs with FEL simulation FAST. Simulation results are compared with experimental data published earlier . Despite experimental data were collected during extended period of time, for the above mentioned publications we selected only those data which corresponded to similar tuning of the accelerator. The same settings of the accelerator were used in the start-to-end simulations.

### 4.1 Radiation energy and fluctuations

Left plot in Fig. 8 presents the average energy in the radiation pulse versus undulator length. Details of experiment are presented in . The interaction length has been changed by means of switching on electromagnetic correctors installed inside the undulator. The value of the orbit kick provided by a corrector was sufficient to stop FEL amplification process downstream the corrector. The radiation energy has been measured by means of an MCP-based detector of 10 mm diameter installed 12 m downstream the undulator . When the FEL interaction is suppressed along the whole undulator length, the detector shows the level of spontaneous emission of about 2.5 nJ collected from the full undulator length. In order to reproduce correctly experimental situation, the simulation results were distorted with the same level of noise (5%) as that provided by the radiation detector. Comparison of experimental and simulation results shows reasonable agreement which is within limits of accuracy of experiment. Experimentally measured power gain length is about 70 cm. Calculations of power gain length for the same electron beam without energy chirp give the value which is almost twice shorter, about 40 cm. Explanation of this puzzle is in strong suppression of the FEL gain by the energy chirp in the electron bunch, about 1% on a scale of cooperation length (to be compared with the FEL parameter $`\rho `$  which is about 0.5%).

Each circle in the left plot in Fig. 8 is the result of averaging over 100 shots. The energy in the radiation pulse fluctuates from shot to shot. The plot for the standard deviation $`\sigma `$, is presented on the right side in this Figure. At the initial stage fluctuations are defined mainly by the fluctuations of the charge in the electron bunch and the accuracy of measurements. When the FEL amplification process takes place, fluctuations of the radiation energy are mainly given by the fundamental statistical fluctuations of the SASE FEL radiation . A sharp drop of the fluctuations in the last part of the undulator is a clear physical confirmation of the saturation process. Detailed measurement of probability distributions were made for the end of the linear regime (undulator length 9 m), and in the nonlinear regime. A comparison of experimental and simulated probability distributions is presented in Fig. 8. It is seen, as expected from theory, that in the linear regime both distributions are gamma-distributions

$$p(E)=\frac{M^M}{\mathrm{\Gamma }(M)}\left(\frac{E}{E}\right)^{M1}\frac{1}{E}\mathrm{exp}\left(M\frac{E}{E}\right),$$

with the value $`M=1/\sigma ^22.6`$ (corresponding to $`\sigma =62`$% at the undulator length 9 m, see Fig:8). Right plot in Fig. 8 shows comparison of the measured and simulated statistical distributions for the nonlinear regime which again are in excellent agreement.

### 4.2 Angular divergence

Figure 9 shows the angular divergence of the FEL radiation. Measurements were performed by means of scanning a 0.5 mm aperture across the radiation beam and recording the transmitted intensity on a downstream detector . The simulated curve has been produced by recalculating the radiation field from the near (undulator exit) into the far zone. The agreement between simulation and measurement is perfect: even fine details in the distribution are well reproduced. We can thus state that our simulation model provides the same field distribution at the undulator exit as it was in the experiment. We should emphase that the disagreement in angular distribution presented in a previous analysis  was a strong indication on rather complicated physical process. Experimental procedure was very reliable, since it was simple relative measurement of the intensity with respect to maximum. Indeed, we can conclude that long spanning tails in the angular distributions is consequence of the strong longitudinal energy chirp which significantly disturbs the beam radiation mode. Experimentally it was not possible to measure spot size of the radiation at the undulator exit. Right plot in Figure 9 shows relevant distribution reconstructed from the simulation data. Figure 9 led us to the conclusion that there should be a high degree of transverse coherence in the FEL radiation, as it was proven in a dedicated experiment at TTF FEL .

### 4.3 Radiation spectra and fluctuations

The next topic of our study is the radiation spectra. Experimentally such spectra can be measured in a single shot by the mean of a monochromator . In such a measurement the procedure is as follows: the photon beam is deflected by plane mirror, passed through a narrow slit, and then dispersed by a grating. The dispersed beam is then detected by CCD camera as illustrated in Fig. 10 in the left-corner plot. Projection of the CCD readings onto $`x`$-axis gives the spectrum (see right-corner plot in Fig. 10). The same procedure was reproduced in the simulations and gave the results presented in the lower plots of Fig. 10 . The range of spectral measurements was rather limited due to limited sensitivity of the spectrometer. It was possible to detect reliably single shot spectra in the nonlinear regime, while in the linear regime only spectra, averaged over many shots are available after procedure of background subtraction. In Figure 12 we present comparison of experimental and simulated spectra for nonlinear regime. We again see not only qualitative, but very good quantitative agreement. Note that the bump on the left slope of the average spectra was a signature for all spectral measurements at TTF FEL starting from the first lasing on February 22, 2000 (see Fig. 12). Such a bump is indeed the consequence of the strong energy chirp along the bunch.

Another important topic of the FEL physics, well documented in the experiment, are the statistical properties of the radiation after narrow-band monochromator. Measurements of fluctuations of the radiation energy after the monochromator have been performed using a narrow-band monochromator of the RAFEL (Regenerative Amplifier FEL ) optical feedback system. The scheme for these measurements is shown in Fig. 15 . The SASE FEL radiation emitted by the electron beam is back-reflected by a plane SiC mirror (RAFEL chamber at the right side of the scheme) onto monochromator (RAFEL chamber at the left side of the scheme). The RAFEL monochromator is a spherical grating in Littrow mounting which disperses the light in the direction of the radiation detector unit (RDU) installed 27 meters downstream. The RDU is equipped with an MCP-based radiation detector with a thin (200 $`\mu `$m) gold wire which plays the role of an exit slit of the monochromator. The design of a spherical grating in Littrow mounting guarantees a resolution of about $`(\mathrm{\Delta }\omega /\omega )_\mathrm{M}10^4`$ which is much less than typical scale of the spike in spectrum (see Fig. 12). Wide dynamic range of MCP-based radiation detector gave the possibility to perform reliable measurements for both the linear and nonlinear regimes.

In the simulation procedure we also took into consideration aperture limitations along the path of the photon beam (9 mm in the undulator, and 6 mm in the spoiler of the electron beam collimator). The reason for this is that diffraction effects at aperture edges lead to mixing of spectrum which initially was strongly correlated with the angle (see Fig. 10).

The plots in Fig. 15 show the probability distributions of the radiation energy after a narrow-band monochromator. SASE FEL operates in the linear regime, the active undulator length during this measurement was 9 m. When SASE FEL operates in linear regime, the probability density must be a negative exponential distribution:

$$p(E)=\frac{1}{E}\mathrm{exp}\left(\frac{E}{E}\right).$$

The solid curves in Fig. 15 represent the negative exponential distribution. So, we can conclude that both, measured and simulated properties well follow the general properties typical for completely chaotic polarized light.

Theory of SASE FEL predicts that in the saturation regime the fluctuations should drop visibly when pulse duration $`T`$ is such that $`\rho \omega T2`$ . At larger values of $`\rho \omega T`$ fluctuations increase and quickly approach 100% level. Since radiation pulse length of TTF FEL is about two cooperation length, we should expect significant suppression of fluctuations in the nonlinear regime, down to 40%. It is seen in Fig. 15 that measured fluctuations drop drastically in the nonlinear mode of operation of TTF FEL. There is not only qualitative agreement, but also quantitative agreement with calculated probability density distribution function. It is worth to mention that such a stabilization of fluctuations is an independent indication for very short pulse duration $`T`$ .

## 5 Discussion

The good agreement between experimental data and simulation results allows us to determine the parameters of the FEL which are not directly accessible experimentally. First of all this refers to the temporal structure of the radiation pulse (see Fig. 16): the computed FWHM pulse duration in saturation regime is about 40 fs, and peak power (averaged over ensemble) is 1.5 GW. We can also conclude that the phase 1 of the SASE FEL at the TTF FEL, was driven by strongly non-Gaussian beam having peak current about 3 kA. Beam dynamics in the accelerator (even at high energy, after bunch compressor) was strongly influenced by the space charge effects.

Our simulations show that rather poor transverse emittance of the electron beam originated from the injector. With a large value of emittance, the beam dynamics is not strongly influenced by CSR effects in the bunch compressor. However CSR effects might be a crucial issue for low-emittance beams as foreseen in future X-ray FELs. Strong longitudinal space charge effect was not expected at TTF since local energy spread before compression (parameter defining spike width and peak current after compression) was strongly overestimated. In future machines this effect can be avoided by an appropriate choice of compression schemes. Therefore, the results of our simulations cannot be directly scaled to more advanced accelerator designs for X-ray FELs.

In this paper we did a ”font-end” comparison: we only compared the results of the simulations with the radiation properties which have been measured reliably and with high accuracy. Although direct comparisons of the beam dynamics simulations with the measurements of the slice parameters of the beam at different location along the accelerator would have been more appropriate, there was impossible due to the lack of adequate beam diagnostics. Our rich experimental experience at TTF shows that tuning of the accelerator for the FEL operation is very difficult task without information about relevant beam parameters. Excellent properties of the FEL radiation were mainly achieved by a global empirical optimization of the machine. Such a procedure might be impossible in much larger and more complicated accelerator systems for X-ray FELs with very large parameter space to be tuned. This underlines the urgent necessity to develop electron beam diagnostics tools which are crucial for the proper operation of X-ray SASE FEL. In particular, a method of the peak current and bunch profile reconstruction, based on detection of infrared coherent radiation from an undulator, would perfectly fit this purpose . Reliable methods for measuring slice emittance and energy spread on a femtosecond time scale should be developed, too.

## Appendix: Parameters of the leading peak in density distribution created in the bunch compressor due to RF nonlinearity

An expression for the beam density distribution after compression, including RF nonlinearity, was derived in . We present here the simplified formulas for the parameters of the leading peak.

Let us consider a long bunch with the constant current $`I_0`$ accelerated off-crest in an RF accelerator with the RF wavelength $`\lambda `$. An energy chirp along the beam is

$$E=E_0\mathrm{cos}(\varphi _0+\mathrm{\Delta }\varphi )$$

Here $`\varphi _0`$ is the phase of a reference particle, and $`\mathrm{\Delta }\varphi =2\pi s/\lambda `$. A position $`s`$ of a particle is positive if it is moving in front of the reference particle. We assume that $`\mathrm{\Delta }\varphi 1`$ and expand the relative energy chirp up to the second order:

$$\stackrel{~}{\delta }=\frac{EE_0}{E_0}\mathrm{\Delta }\varphi \mathrm{tan}\varphi _0\frac{(\mathrm{\Delta }\varphi )^2}{2}$$
(1)
For a Gaussian uncorrelated energy spread a distribution function in the vicinity of the reference particle can be described as follows:

$$f(s,\delta )=\frac{I_0}{\sqrt{2\pi }\sigma _\delta }\mathrm{exp}\left[\frac{\left(\delta +2\pi s\mathrm{tan}(\varphi _0)/\lambda +2\pi ^2s^2/\lambda ^2\right)^2}{2\sigma _\delta ^2}\right],$$
(2)
where $`\sigma _\delta `$ is the relative rms energy spread. The normalization is chosen such that after integration over $`\delta `$ we get the current.

Behind the bunch compressor a particle position $`s_f`$ (with respect to a nominal particle) is connected with its position before compression $`s_i`$ and energy deviation $`\delta `$ as

$$s_fs_i+R_{56}\delta +T_{566}\delta ^2+\mathrm{},$$
(3)
where $`R_{56}`$ and $`T_{566}`$ are the first and the second order momentum compaction. In our consideration we choose the signs such that for a bunch compressor chicane $`R_{56}>0`$ and $`T_{566}3R_{56}/2`$ (for a small bending angle). It is seen from (1) and (3) that the phase space distribution after compression has a parabolic shape. The reference particle is positioned at the fold-over when

$$\mathrm{tan}(\varphi _0)=\frac{\lambda }{2\pi R_{56}}$$

Also, analyzing Eqs. (1) and (3), we come to the conclusion that the contribution of $`T_{566}`$ term can be neglected under the condition $`\varphi _01`$. We also assume that $`\sigma _\delta \varphi _0^2`$. The final distribution function is obtained from (2) by a simple substitution $`ssR_{56}\delta `$ and is simplified with the help of the above mentioned assumptions:

$$f(s,\delta )=\frac{I_0}{\sqrt{2\pi }\sigma _\delta }\mathrm{exp}\left[\frac{\left(s/R_{56}+2\pi ^2R_{56}^2\delta ^2/\lambda ^2\right)^2}{2\sigma _\delta ^2}\right].$$
(4)
To obtain the dependence of current on s after compression one should make the integration:

$$I(s)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}\delta f(s,\delta )$$

After normalization we get:

$$I(\widehat{s})=I_0Cg_1(\widehat{s}),$$
(5)
where

$$\widehat{s}=\frac{s}{R_{56}\sigma _\delta },C=\frac{\lambda }{\pi R_{56}\sqrt{2\sigma _\delta }},$$

and the function $`g_1(\widehat{s})`$ is

$$g_1(\widehat{s})=\sqrt{\frac{2}{\pi }}\underset{0}{\overset{\mathrm{}}{}}x\mathrm{exp}\left[\frac{(\widehat{s}+x^2)^2}{2}\right].$$
(6)
The plot of the function $`g_1(\widehat{s})`$ is presented in Fig. 17. The maximal value of the function $`g_1(\widehat{s})`$ is close to 1, $`\mathrm{max}(g_1(\widehat{s}))1.02`$, so that $`C`$ describes the enhancement of current. The full width (at half maximum) of this curve is 4.8, or in dimensional notations:

$$(\mathrm{\Delta }s)|_{\mathrm{FWHM}}4.8R_{56}\sigma _\delta $$
(7)
At the level of 0.8 the full width $`\mathrm{\Delta }\widehat{s}`$ is equal to 2. We can also estimate the length of the beam slice before compression, contributing to the leading peak after compression \- it is of the order of $`\lambda \sqrt{\sigma _\delta }`$. Thus, if the current (before compression) only weakly changes on this scale, our approximation of constant current is valid.

It is interesting to know the slice energy spread after compression. It is easy to obtain from (4) the following result:

$$\delta _{\mathrm{rms}}(\widehat{s})=\sigma _\delta Cg_2(\widehat{s}),$$
(8)
where

$$g_2(\widehat{s})=\left[\frac{\underset{0}{\overset{\mathrm{}}{}}xx^2\mathrm{exp}\left[(\widehat{s}+x^2)^2/2\right]}{\underset{0}{\overset{\mathrm{}}{}}x\mathrm{exp}\left[(\widehat{s}+x^2)^2/2\right]}\right]^{1/2}$$
(9)
The plot of this function is presented in Fig. 18 In the slice with maximum current, $`\widehat{s}=0.8`$, the function $`g_2`$ takes the value of 0.89. The ratio111We do not call this quantity ”longitudinal brightness” to avoid a possible confusion since the ratio of local current to slice energy spread does not have to be conserved. Instead, the phase space density is conserved. Note also that in our simplified treatment we did not take care of the simplecticity of transformation. Nevertheless, this practically does not influence our results.

$$\frac{I(s)}{\delta _{\mathrm{rms}}(s)}=\frac{I_0}{\sigma _\delta }\frac{g_1(\widehat{s})}{g_2(\widehat{s})}$$

can be important for FEL operation. In Fig. 19 we plot the function $`g_1(\widehat{s})/g_2(\widehat{s})`$. One can see some enhancement of the ratio of local current to slice energy spread with respect to uncompressed beam. The full width of the curve is about 3, beyond this good range it quickly goes down, so that only a small part of the beam is favorable for lasing (if it is not spoiled by collective effects). Note that rms energy spread is an adequate quantity for SASE FEL calculation when it is much less than the FEL parameter $`\rho `$ . When it becomes comparable to $`\rho `$, one should take into account the actual distribution (to know only dispersion is not sufficient). When it is much larger that $`\rho `$, one can think of two independent beams, having different energies and twice lower current.

Finally, let us present a numerical example. An electron beam with the current 80 A and rms energy spread 3.5 keV is accelerated up to 135 MeV in the RF accelerator with $`\lambda =23`$ cm. After that it is compressed in a bunch compressor with $`R_{56}=23`$ cm. The peak current after compression is then 3.6 kA, rms energy spread in the slice with maximum peak current is about 140 keV, and full width of the leading peak is 28 $`\mu `$m.