text:
# Triplet correlations in two-dimensional colloidal model liquids

## 1 Introduction

The static structure of a simple fluid is commonly described in terms of the $`m`$-body distribution functions $`g^{(m)}`$, measuring the probability densities of finding two, three, and more particles at specified positions in space. The present paper is concerned with three-body distribution functions, as obtained from measured configurations in two-dimensional colloidal model fluids.

There are many examples in statistical mechanics where triplet correlations in classical fluids are important. They play an essential role, for example, in earlier theories of the pair correlation function of a fluid, such as the Born-Green-Yvon integral equation theory, which explicitly shows the relationsship between pair and triplet correlation functions . Triplet correlations are furthermore used in perturbation theories for static fluid properties , in theories of transport properties , but also to describe solvent reorganization processes around solutes  or systems under shear-flow . And, finally, one is again led to triplet correlation functions if one is interested in the temperature or density derivative of the pair-correlation function $`g^{(2)}(r)g(r)`$ .

For these and other reasons, it is not surprising to see how often triplet correlations have been theoretically investigated, in the overwhelming number of cases, by computer simulation studies. Of central importance in these studies is Kirkwood’s idea  to approximate the triplet distribution function by three different pair distribution functions (the Kirkwood superposition approximation (KSA) explained in more detail below). The reliability and quality of the KSA has been tested extensively for pure hard-sphere systems , for Lennard-Jones (LJ) fluids , for liquid sodium , but also for Yukawa fluids , with the general result that all qualitative features of $`g^{(3)}`$ are usually well described by the KSA. However, quantitative differences at small distances and high densities can be appreciable, see, for instance, . Mc Neil et al.  compared their simulations also with analytical theories going beyond the KSA. More recent simulations aiming at an understanding of higher-order correlations are performed in 2-component asymmetric electrolyte , in an aqueous 1:1 electrolyte , and in strong 3:3 and 1:3 electrolytes . Finally, one has to mention theoretical work done using integral equation theory , see  for a review.

All these theoretical efforts contrast with the situation on the experimental side where there are only very few papers dealing with three-body correlations. One indirect way to obtain information on $`g^{(3)}`$ is via the isothermal pressure derivative of the fluid structure factor $`S(q)/P`$ which is related to the triplet distribution function . This idea has been exploited in rare-gas fluids mainly by Egelstaff and coworkers in a series of experimental papers: in Krypton , Argon , Helium  and Neon . But, to our knowledge, there is no direct measurement of $`g^{(3)}`$. It is clear where the problems come from: x-ray or neutron diffraction data provide information just on the fluid structure factor, i.e., essentially on the pair-correlation function $`g(r)`$. A measurement of $`g^{(3)}`$, however, requires the knowledge of the positions of three particles at the same time which is technically very demanding to obtain in 3D samples (see however ).

Here lies the advantage of video-microscopy, a modern experimental technique applied to colloidal systems to directly measure all particles’ positions at all times. Essentially, one follows the phase-space trajectory of the system as it evolves in time, and has thus the same amount of information as one obtains in a simulation. It is rather straightforward to calculate higher-order correlation functions from the measured configurations. In , we have recently realized this idea and calculated three-body correlation functions from configurations obtained by video-microscopy from a two-dimensional system of magnetic colloids. This is probably the first direct measurement of a three-body correlation function. In the present paper now we give details not mentioned in , outline a few technical points and present calculations of $`g^{(3)}`$ also for 2D Yukawa systems.

## 2 Theoretical background

### 2.1 Distribution functions

The radial pair and three-body distribution functions $`g(r)`$ and $`g^{(3)}(𝐫_1,𝐫_2,𝐫_3)`$ are used to describe the local structure of a fluid and can best be defined by means of the molecular distribution functions $`n^{(m)}(𝐫^m)`$ . The lowest order member of this class of functions, $`n^{(1)}(𝐫_1)dr_1`$, is the probability of finding a particle in the volume $`dr_1`$ at $`𝐫_1`$. In an homogeneous medium it is equal to the density $`\rho `$. $`n^{(2)}(𝐫_1,𝐫_2)dr_1dr_2`$ is the joint probability of finding one particle in the volume $`dr_1`$ at $`𝐫_1`$ and a second in the volume $`dr_2`$ at $`𝐫_2`$ and so on for higher orders. Dividing out the asymptotic dependency we obtain the set of $`m`$-body distribution functions $`g^{(m)}(𝐫^m)`$:

$$g^{(m)}(𝐫^m)=\rho ^mn^{(m)}(𝐫^m).$$
(1)
The first non-trivial radial distribution function is $`g(𝐫_1,𝐫_2)`$ which in a homogeneous, isotropic system depends only on the particle separation $`r=|𝐫_2𝐫_1|`$. Next in order is $`g^{(3)}(𝐫_1,𝐫_2,𝐫_3)`$ describing the probability of finding triplets of particles. Configurations of three particles are uniquely characterized by three independent parameters (in a homogeneous, isotropic medium), which can be chosen to be the distances between the particles $`r=|𝐫_2𝐫_1|`$, $`s=|𝐫_3𝐫_1|`$ and $`t=|𝐫_3𝐫_2|`$. Therefore, $`g^{(3)}=g^{(3)}(r,s,t)`$.

The potential of mean force is defined by

$$\beta w^{(m)}=\mathrm{ln}g^{(m)}$$
(2)
with $`\beta =1/kT`$ for the inverse thermal energy. The potentials of mean force have to be distinguished from the direct potentials; so, $`w^{(2)}(r)`$ should not be confused with the direct pair-potential $`u(r)`$ in the system. The difference between both quantities can best be understood by means of the Born-Green equation ,

$$\frac{w^{(2)}(r_{12})}{𝐫_1}\frac{u(r_{12})}{𝐫_1}=\rho \frac{u(r_{13})}{𝐫_1}\frac{g^{(3)}(𝐫_1,𝐫_2,𝐫_3)}{g(r_{12})}𝑑𝐫_3,$$
(3)
relating the difference between the mean force and the direct pair-force to an integral over the force on particle 1 at $`𝐫_1`$ due to a particle at $`𝐫_3`$, weighted by the probability $`\rho g^{(3)}d𝐫_3/g(r_{12})`$ of finding a particle in $`d𝐫_3`$ at $`𝐫_3`$ when it is known that other particles are located at $`𝐫_1`$ and $`𝐫_2`$. This equation is a member of the Born-Green-Yvon (BGY) hierarchy and is exact if pairwise interactions can be assumed. Inserting $`\beta w^{(2)}=\mathrm{ln}g`$ and using the KSA as a closure relation, eq. (3) yields the BGY integral equation for $`g(r)`$ .

### 2.2 Kirkwood Superposition Approximation

In a dilute system the overall interaction is dominated by the interaction of individual pairs. Therefore the probability of finding a certain arrangement of particles is nothing but the joint probability of finding individual pairs. For a triplet with particle distances $`r`$, $`s`$ and $`t`$ the approximate triplet distribution function $`g_{SA}^{(3)}(r,s,t)`$ is then the product of the pair distribution functions

$$g_{SA}^{(3)}(r,s,t)=g(r)g(s)g(t).$$
(4)
This is the so called Kirkwood superposition approximation (KSA) . Introducing the factor $`G`$ correcting the error made by the KSA

$$g^{(3)}(r,s,t)=g_{SA}^{(3)}(r,s,t)G(r,s,t)$$
(5)
and the correction potential of mean force $`\mathrm{\Delta }w^{(3)}`$

$$\mathrm{\Delta }w^{(3)}(r,s,t)=\mathrm{ln}G(r,s,t)/\beta $$
(6)
we obtain for the triplet potential of mean force $`w^{(3)}`$

$$w^{(3)}(r,s,t)=w^{(2)}(r)+w^{(2)}(s)+w^{(2)}(t)+\mathrm{\Delta }w^{(3)}(r,s,t).$$
(7)
All pair correlations in $`g^{(3)}`$ are included in $`g_{SA}^{(3)}`$ while $`G`$ quantifies the extent of intrinsic correlations due to the simultaneous presence of a triplet of particles. For this reason, $`G`$ is called the triplet correlation function. In other words, if $`G`$ is unity everywhere, there are no genuine triplet correlations in the system, but only the trivial ones that can be expressed by pair correlation functions. $`\mathrm{\Delta }w^{(3)}`$ therefore measures the extra correlation energy of three correlated particles relative to the energy of superposed correlated pairs. One should be careful not to confuse $`\mathrm{\Delta }w^{(3)}`$ with a real three-body potential $`u^{(3)}`$: $`\mathrm{\Delta }w^{(3)}`$ is a correlation energy and can thus have non-vanishing values even in system with particles interacting exclusively via pair-wise additive potentials $`u(r)`$, i.e., in cases where higher-order potentials such as $`u^{(3)}`$ are strictly zero.

## 3 Experimental systems and technical remarks

### 3.1 Two colloidal model systems

The triplet functions that we here present are based on a large number of particle configurations obtained from two different colloidal model systems by means of video-microscopy. Both systems are two-dimensional (or more precisely: quasi-2D). In the first system particles interact via a tuneable magnetic dipole-dipole interaction and in the second via an electrostatic double-layer interaction.

The first system consists of paramagnetic spherical colloidal particles with an diameter of $`4.7\mu m`$. They are located at the bottom of an hanging water droplet whose surface can be accurately controlled to be almost perfectly flat. The positions of the particles are recorded using digital video-microscopy with subsequent image-processing on the computer. The field of view has a size of $`520\times 440\mu m`$ containing typically about $`10^3`$ particles. A magnetic field $`B`$ applied perpendicular to the air/water interface induces in each particle a magnetic moment $`M=\chi B`$ which leads to a repulsive dipole-dipole pair-interaction energy of

$$\beta u(r)=\frac{\mathrm{\Gamma }}{(\sqrt{\pi \rho }r)^3}$$
(8)
with the interaction strength given by $`\mathrm{\Gamma }=\beta (\mu _0/4\pi )(\chi B)^2(\pi \rho )^{3/2}`$. This is the only relevant contribution to the interparticle-potential which is hence conveniently and reversibly adjustable by varying $`\mathrm{\Gamma }`$ through the external field $`B`$. Note in particular that by introducing $`\mathrm{\Gamma }`$ we have scaled out the density. $`\mathrm{\Gamma }`$ is thus the only parameter determining the phase-behavior of the system: for $`\mathrm{\Gamma }<57`$ the system is liquid, for $`\mathrm{\Gamma }>60`$ it is solid, and in between, i.e. for $`57<\mathrm{\Gamma }<60`$, it shows an hexatic phase . The system can be regarded as an almost ideal 2D model system as the out-of-plane motion of the particles corresponds to less than 1 % of their diameter. Details about the preparation of the samples and the experimental set-up can be found in .

The second system is an aqueous suspension of highly charged sulphate-terminated PS particles of $`\sigma =3\mu m`$ diameter, confined between two glass plates with a 1 mm spacing. The particles were furthermore exposed to vertical light forces which pushed them toward the negatively-charged silica plate at the bottom of the cell, confining the system more effectively to two dimensions. The particle center positions were analyzed on-line with an imaging processing software. The particle density was varied in the system by a scanned optical laser tweezer which acts as a corral for the investigated 2D colloidal suspension. The particles interact with a screened Coulomb potential

$$\beta u(r)=\left(\frac{Z}{1+\kappa \sigma /2}\right)^2e^{\kappa \sigma }\lambda _B\frac{e^{\kappa r}}{r},$$
(9)
with the (effective) particle charge $`Z`$, the particle diameter $`\sigma `$, the Bjerrum-length $`\lambda _B`$ and the inverse Debye-screening-length $`\kappa `$. Via the inversion of the Ornstein-Zernicke equation for a low density system we can determine the potential quite accurately. This experiment and also the inversion procedure has been described in great detail in . In fact, we here compute triplet correlation functions based on exactly the same set of configurations as already evaluated for other purposes in . Pair-correlation functions corresponding to the triplet correlation functions of this work can also be found in .

In both systems, we made sure that the system is well equilibrated and used about 200 statistically independent configurations with approximately 500-1000 particles. We also performed standard Monte-Carlo (MC) simulations, which offers us a third way to generate the necessary set of particle configurations. The simulations produced typically 500 configurations with 2000 particles yielding much better statistics and the opportunity to check the experimental results for statistical errors.

### 3.2 How to compute triplet distribution functions

The numerical procedure for computing $`g^{(3)}`$ is essentially the same as the one used in simulations to determine $`g(r)`$. For every configuration we counted the triplets with sidelengths $`r,s`$ and $`t`$ in an array $`[R_i,S_j,T_k]`$ such that $`R_i\mathrm{\Delta }R/2r<R_i+\mathrm{\Delta }R/2`$ and accordingly for $`s`$ and $`t`$ with $`\mathrm{\Delta }R=\mathrm{\Delta }S=\mathrm{\Delta }T`$. While normalizing $`g`$ means each slot is divided by the area $`A(R_i)=2\pi R_i\mathrm{\Delta }R`$ it covers, we can write down in analogy an analytic expression for $`g^{(3)}`$ valid for infinitesimal small $`dR`$:

$$A_{an}(R_i,S_j,T_k)=\frac{4\pi R_iS_jT_kdR^3}{\sqrt{R_i^2S_j^21/4(R_i^2+S_j^2T_k^2)^2}}.$$
(10)
This expression also works fine for finite $`\mathrm{\Delta }R`$ (replace $`dR`$ by $`\mathrm{\Delta }R`$ in eq. (10)) except in the case of triangles with one angle close to 180 degrees. For example, a triangle with $`r=1.4\mathrm{\Delta }R,s=1.4\mathrm{\Delta }R,t=2.6\mathrm{\Delta }R`$ would be sorted into the slots $`i=1,j=1,k=3`$. The center values of these slots obviously do not represent a triangle and thus the normalizing procedure fails because the square root in eq. (10) turns imaginary. While we can neglect these slots for most of the plots presented, we cannot ignore these counts when integrating over the whole $`g^{(3)}`$ as it is done in the section on the Born-Green equation. In these cases we are not allowed to replace $`dR`$ by $`\mathrm{\Delta }R`$ in eq. (10), but, instead, we have to compute the normalization factor numerically. For a fixed distance $`r`$ one can easily calculate the intersection area $`A(r,S_j,T_k)`$ of two circular rings with radii $`S_j,T_k`$ and width $`\mathrm{\Delta }R`$. The numerical normalization factor $`A_{num}(R_i,S_j,T_k)`$ is then

$$A_{num}(R_i,S_j,T_k)=2\pi \underset{R_i\mathrm{\Delta }R/2}{\overset{R_i+\mathrm{\Delta }R/2}{}}rA(r,S_j,T_k)𝑑r.$$
(11)
More technical details can be found in the appendix of .

## 4 Results: Three-body distribution functions

We first concentrate on the equilateral triangle configuration. For this special configuration we have $`r=s=t`$ and the superposition approximation reduces to $`g_{SA}^{(3)}=g(r)g(r)g(r)=(g(r))^3`$. Therefore, $`(g_{SA}^{(3)}(r,r,r))^{1/3}=g(r)`$ and a comparison of $`g(r)`$ and $`(g^{(3)}(r,r,r))^{1/3}`$ reveals the importance of the triplet correlation function $`G(r,r,r)`$ in eq. (5). Fig. (1) shows $`(g_{SA}^{(3)}(r,r,r))^{1/3}`$ and $`(g^{(3)}(r,r,r))^{1/3}`$ for the system with the magnetic colloids (Fig. (1.a)) and that with the charged colloids (Fig. (1.b)). For each system type, we have analyzed three different densities ($`\mathrm{\Gamma }`$-values), where the systems are always deep in the liquid phase. These are: $`\mathrm{\Gamma }=4,14,46`$ for the magnetic system and $`\rho \sigma ^2=0.037,0.167,0.186`$ for the charged system.

As expected, at low density (small $`\mathrm{\Gamma }`$) Kirkwood’s superposition approximation holds very well. These are then weakly interacting systems in which particles sparsely ”meet” other particles and usually interact only with a single other particle. The probability of finding a certain triplet is thus nothing but the joint probability of finding the appropriate pairs, and therefore the KSA is justified by the independence of the particular pair probabilities. Increasing the interaction strength $`\mathrm{\Gamma }`$ (i.e., the density in the charged system) particles begin to interact with more than one other particle at a time. The assumption that a third particle will not interfere with the statistical distribution of the second particle begins to deteriorate. One can see that in Fig. (1): the first peak of $`g^{(3)}(r,r,r)^{1/3}`$ is always reduced in height compared to $`g(r)`$ and the first minimum deviates even stronger. It is always more pronounced for $`g^{(3)}(r,r,r)^{1/3}`$ and it changes its shape from the sine-like oscillation of $`g(r)`$ to a more asymmetric one. This seem to be just a packing effect which is relatively independent of the interaction potential between the particles, as is evident from the surprisingly close resemblance between the curves in Fig. (1.a) and (1.b). Note, however, that shape and height of the first peak is different for both systems, as one would expect for different pair-potentials. The dashed and solid lines represent the results of our simulations, using the pair-potentials in eq. (8) and (9). In the charged system, the prefactor and the inverse screening length has been used as a fit-parameter. The good agreement between simulation and experimental data found for both the pair- and triplet distribution functions demonstrates that the correlations can in both systems be understood in terms of a picture of pair-wise interacting particles, although it has to be emphasized that in the charged system the pair-potential is an effective (density-dependent) one in which – through the fitting procedure – many-body contributions are absorbed into the pair-potential, see discussions in . By contrast, many-body interactions in the magnetic system can be safely ruled out .

Another way to visualize $`g^{(3)}(r,s,t)`$ is to vary two variables while fixing one variable which we here choose to be the variable $`r`$. We take $`r=r_{max}^{(2)}`$ which is the distance where $`g(r)`$ has its first peak (almost equal to $`1/\sqrt{\rho }`$). It is the most probable distance between any pair of particles in the fluid. For this choice the KSA leads to

$$g_{SA}^{(3)}(r_{max}^{(2)},s,t)=g(r_{max}^{(2)})g(s)g(t)$$
(12)
and for $`s\mathrm{}`$ and $`t\mathrm{}`$ $`g_{SA}^{(3)}`$ goes to the maximum value of pair correlation function, $`g_{SA}^{(3)}(r_{max}^{(2)},s,t)g(r_{max}^{(2)})`$. To be able to plot $`g^{(3)}`$ directly in the $`(x,y)`$-plane we have to transform the variables $`s`$ and $`t`$ accordingly, $`s=s(x,y)`$ and $`t=t(x,y)`$. Fig. (2) shows $`g_{SA}^{(3)}(r_{max}^{(2)},s,t)`$ for the $`\mathrm{\Gamma }=46`$ measurement of the magnetic system in the $`(x,y)`$ plane (co-ordinate system as defined in the plot). The lighter grey to white areas correspond to values below and the darker grey to black areas to values above the limiting value $`g(r_{max}^{(2)})`$. Clearly, the plot is highly redundant, for symmetry reasons it would have been sufficient to just show one quarter of it. In addition, we show circles with radii corresponding to the first and second maximum (white circles) and the first minimum (black circles) of $`g(r)`$ around each particle.

With these circles the KSA prediction can be easily understood. First we expect a very pronounced peak where the two circles with the shortest radii intersect. At this position is $`s=r_{max}^{(2)}`$ and $`t=r_{max}^{(2)}`$; therefore both $`g(s)`$ and $`g(t)`$ are at their maximum and that maximizes $`g_{SA}^{(3)}(r_{max}^{(2)},s,t)`$ in eq. (12). This position corresponds to an equilateral triangle configuration. Accordingly, the first minimum on the $`y`$-axis occurs where the two black circles intersect. Following the inner white circle around the right particle means keeping $`s=r_{max}^{(2)}`$ constant. At the point where this line intersects the black circle – indicating the first minimum of $`g(t)`$ – we encounter a local minimum along that line. Following further the circle $`s=r_{max}^{(2)}`$ the second maximum of $`g(t)`$ becomes important. But since the two circles intersect almost tangential there is no distinctive structure but a plateau-like ”banana” enclosing the particle on the outside. In a similar way, all other local maxima and minima in $`g_{SA}^{(3)}`$ can be explained by combining the maxima and minima of $`g(t)`$ and $`g(s)`$.

Fig. (3) compares the full triplet correlation function $`g^{(3)}`$ (right half of each figure) with the triplet correlation function in the KSA, $`g_{SA}^{(3)}`$ (left half of each figure). Again, we consider both the magnetic and the charged system, at the highest density (highest $`\mathrm{\Gamma }`$). To show the structure of $`g^{(3)}`$ more clearly we have set the lower end of the grey-level code to the limiting value $`g(r_{max}^{(2)})`$, so only the peaks, i.e. values where $`g^{(3)}>g(r_{max}^{(2)})`$, are shown. The stripes that can be seen especially close to the x-axis result from the transformation $`g^{(3)}(r_{max},s,t)`$ to $`g^{(3)}(r_{max},x,y)`$ and appear due to limited statistics. An hexagonal lattice with a lattice constant $`a=r_{max}^{(2)}`$ is superposed. As in all our plots, we plot the cubic root of $`g^{(3)}`$, to make the resulting numbers comparable to $`g(r)`$.

First of all one notices that the ”banana”-like structure is not plateau-like anymore but shows a separation into three distinctive peaks. For $`\varphi =120^{}`$ (the angle $`\varphi `$ is defined in Fig. (3)) the peak matches almost exactly the grid point of the underlying hexagonal lattice while for $`\varphi =180^{}`$ is slightly shifted to the inside and is not so pronounced. In the shell of second nearest neighbors there is also some inner structure developing but a clear correspondence to the lattice points is not yet recognizable. By contrast, there is no such correspondence between the peak-structure of $`g_{SA}^{(3)}`$ and the lattice sites when the KSA is applied. So this crystal-like local order that shows up already between particles in the fluid phase, is obviously a genuine effect of triplet correlations. Again, there is a remarkably close resemblance between the results for the magnetic and the charged system, demonstrating that these correlation effects have little to do with the properties of the pair-potential.

To quantify our analysis, we now concentrate on the first neighbor shell around the two central particles in Fig. (3), and fix $`r=r_{max}^{(2)}`$ and $`s=r_{max}^{(2)}`$ while changing $`t=t(\varphi )`$ by varying the angle $`\varphi `$ (defined in Fig. (3)) between $`0^{}`$ and $`180^{}`$. Fig. (4) shows $`g_{SA}^{(3)}(r_{max}^{(2)},r_{max}^{(2)},t(\varphi ))`$ and $`g^{(3)}(r_{max}^{(2)},r_{max}^{(2)},t(\varphi ))`$ as a function of $`\varphi `$, for the magnetic system. The corresponding results for the charged system look again almost identical, and need not be further considered. Varying $`\varphi `$ and thus $`t(\varphi )`$ while fixing $`r=r_{max}^{(2)}`$ and $`s=r_{max}^{(2)}`$, we pass through all lattice points of the particle’s first coordination shell (arrows in Fig. (4) mark positions of lattice points). It can be clearly seen that at up to $`\mathrm{\Gamma }=14`$ $`g^{(3)}`$ and $`g_{SA}^{(3)}`$ are almost identical and that at $`\mathrm{\Gamma }=46`$ the full correlation function $`g^{(3)}`$ shows a peak-structure in close correspondence with the hexagonal structure. This is not the case for the the approximated triplet distribution function $`g_{SA}^{(3)}`$.

Also given in Fig. (4) is the function $`\mathrm{ln}G`$, i.e. $`\mathrm{\Delta }w^{(3)}`$, of eq. (6). It is evident how $`\mathrm{\Delta }w^{(3)}`$ gradually forms on increasing $`\mathrm{\Gamma }`$, with values up to one $`kT`$ . It is also seen that the regions of attractive and repulsive correlation energies $`\mathrm{\Delta }w^{(3)}`$ correspond to the correcting effect which the function $`G`$ has on $`g_{SA}^{(3)}`$ to ensure that $`g^{(3)}`$ adapts locally to the hexagonal symmetry. We conclude that it is an effect entirely due to three-particle correlations, i.e. due to the function $`G`$, which is responsible for the observed formation of a crystal-like local environment around particles in a fluid well below the freezing transition. All quantities displayed in Fig. (4) are also compared to the corresponding results of our MC simulations, showing again a remarkably good agreement.

While Fig. (4) follows $`g^{(3)}`$ along a circle around the right particle in Fig. (3), Fig. (5) represents a cut in the $`y`$-direction ($`x`$ = 0) in Fig. (3). In this direction we have generally found the highest triplet correlation energies $`\mathrm{\Delta }w^{(3)}`$. For $`r=r_{max}^{(2)}`$ $`\mathrm{\Delta }w^{(3)}`$ can become as high as 2 $`kT`$, and if $`r`$ is slightly increased to $`1.33r_{max}^{(2)}`$, we find an even higher correlation energy of more than 4 $`kT`$, see Fig. (5.b). Fig. (6) shows a sequence of pictures similar to that in Fig. (3). While in Fig. (3) we fixed $`r`$ to $`r_{max}^{(2)}`$, $`r`$ is now varied. We expect triplet correlations to vanish if all three distances $`r`$, $`s`$ and $`t`$ become large, and that then $`g^{(3)}g_{SA}^{(3)}`$. Indeed, one observes from Fig. (6) that both halves of each figure (displaying $`g^{(3)}`$ and $`g_{SA}^{(3)}`$) show more and more resemblance the larger $`s`$ and $`t`$ are, i.e., in a region far away from the central pair, and that the extent of this region increases with increasing $`r`$. For example, for $`r=2r_{max}^{(2)}`$, differences between both sides are to be seen only on the midplane.

The closer the system approaches the crystalline phase the more pronounced do we expect the lattice points to be occupied. This is tested in Fig. (7) which compares the triplet distribution function for the system in the solid phase ($`\mathrm{\Gamma }=80`$) with that for the system in the liquid phase ($`\mathrm{\Gamma }=46`$). For the solid phase, we performed MC simulations, starting from a perfect hexagonal lattice, while the $`\mathrm{\Gamma }=46`$ distribution is the same as in the right half of Fig. (3.a). With regard to the correlations between the central pair and the first coordination shell, there is hardly any difference between the liquid and the solid phase. Pronounced differences are observable, however, in the second shell: in the liquid phase the next nearest neighbors are broadly distributed midway between adjacent lattice nodes, while the $`\mathrm{\Gamma }=80`$ distribution clearly correlates much better with the lattice structure (see, for example, the lattice point in the second shell next to the mid-plane between the central pair). However, even for $`\mathrm{\Gamma }=80`$ this correspondence is far from perfect; especially close to the $`\varphi =180^o`$ direction there is still an extended smeared-out distribution showing no clear preference for certain lattice points. Clearly, approaching $`T0`$ ($`\mathrm{\Gamma }\mathrm{}`$), one will ultimately observe peaks in $`g^{(3)}`$ positioned exclusively on the lattice points.

Krumhansl and Wang used eq. (3) to check the accuracy of the KSA in their early simulations of a LJ fluid . We copy their idea, but use experimental data to demonstrate the importance of three-particle correlations. We numerically computed the right-hand side of eq. (3) using both the full and the approximated triplet function, $`g^{(3)}`$ and $`g_{SA}^{(3)}`$, of the $`\mathrm{\Gamma }=4`$ and $`\mathrm{\Gamma }=46`$ measurement and compare it in Fig. (8) to the left-hand side of eq. (3), evaluated using $`u(r)`$ and $`g(r)`$. For the weakly interacting system ($`\mathrm{\Gamma }=4`$), $`G`$ is unity in value everywhere: triplet correlations are unimportant, and, accordingly, the Born-Green equation can be satisfied with triplet functions based on the KSA, see Fig. (8). In the strongly interacting system ($`\mathrm{\Gamma }=46`$), however, the KSA fails completely. Three-particle correlations have to be taken fully into account to obtain the correct difference between mean and direct force via the Born-Green equation.

A number of other approximations for $`g^{(3)}`$ are known , and could be checked in a similar way. This includes the Schofield equation , relating $`g(r)/\rho `$ to $`g^{(3)}`$, which in turn is the basic equation for a number of thermodynamic consistency relationships . From our results, we can expect that, here again, $`g^{(3)}`$ at high $`\mathrm{\Gamma }`$ cannot be approximated, but has to be taken in its full form. Our results confirm predictions from simulation studies on LJ and hard-sphere systems showing that qualitatively the structure of $`g^{(3)}`$ is correctly described by the KSA, but that quantitative failures can be appreciable .

## 5 Conclusion and closing remarks

We used experimental data obtained from digital video microscopy of particles in a two-dimensional colloidal model fluid. From the measured configurations, we subsequently calculated particle correlation functions; in particular, triplet correlation functions, containing much more information on the relative spatial arrangement of particles than the pair-correlation function. We studied two colloidal systems: one with particles interacting via a Yukawa pair-potential, and another with particles interacting with a $`\mathrm{\Gamma }/r^3`$ potential. We have found very similar results for both systems: When the liquid is near the freezing transition (high density, high $`\mathrm{\Gamma }`$), the deviations of the three-particle correlation function from unity are considerable (with correlation energies as high as 4 $`kT`$). We observe the formation of a crystal-like local environment around particles in a fluid well below the freezing transition. Clearly, these are packing effects as is evident from the close resemblance between the results for the two systems having completely different pair-interaction potentials. All quantities examined have also been compared to results from MC simulations, and showed always good agreement. We have finally used our experimental data in combination with the Born-Green equation to demonstrate that triplet correlations are also significant when integrated over the whole volume. For the strongly interacting magnetic system ($`\mathrm{\Gamma }=46`$), the KSA has been shown to fail completely. Three-particle correlations are thus seen to be important not only to obtain locally the correct structure, but also to obtain globally the correct difference between mean and direct force via the Born-Green equation.

In principle, higher-order correlation functions could be computed from the measured configurations in very much the same way, but it is not clear to us what they could be useful for. We should also remark that it soon will be possible to analyze higher-order correlation functions also in 3D colloidal samples, with positional data recorded using modern confocal microscope techniques. One has to be aware, however, that in three dimensions a similar correspondence between the peaks in $`g^{(3)}`$ and an underlying crystal lattice should be much harder to find. In 3D, every triplet of particle lies, of course, also in a plane, and can accordingly be plotted as in Fig. (2) and (4). However, then there is not one, but a superposition of many possible lattice planes that one has to compare this distribution with. In this respect, the role of a crystal lattice in determining the structure of a liquid is probably more pronounced in 2D than it is in 3D.

HHvG and CR wish to thank Matthias Brunner and Prof. Clemens Bechinger for all the many intensive discussions on this and other topics, for their help and assistance, and the steady and close cooperation. We gratefully acknowledge continuous support and stimulating discussion with Prof. Rudolf Klein, and financial support from the Deutsche Forschungsgesellschaft through the Konstanz SFB 513.

## References


format: text

text:
# Analysis of laboratory nucleosynthesis products

## Introduction

This work summarizes the results of analysis of the nuclear transmutation products in the supercompressed substance obtained in 1999–2003 by the Electrodynamics Laboratory ‘Proton–21’ (Ukraine, Kiev) that was established to develop a novel approach for radioactive waste utilization.

The project was based on the idea that it is possible to create a special kind of disturbance impact on a solid substance that will induce, under certain conditions, the self–organizing avalanche process of supercompression up to the collapse state at which the conditions for collective multiparticle transmutation reactions of radioactive elements into stable ones will become possible, owing to the redistribution of neutrons in the substance volume and increase of their concentration at the center of the collapse. A few experimental setups have been built in the laboratory. The first successful experiment was performed on February 24, 2000. By March, 2003, 4037 dynamic compressions of solid targets were performed leading to specific explosions and the radial dispersion of a transformed target material from the collapse zone.

## I Multinuclear reactions in a superdense substance

### I.1 Experiment

The main experimental setup realizes the consecutive space–time energy compression and is capable to impact on the solid target up to 1 kJ of energy during the pulse time less than 100 ns. The electron beam is used as an initial carrier of the concentrated energy. On the final stage, the substance of target is compressed up to a density of $`>10^{26}`$ cm-3, while the estimated power density in the collapse volume exceeds $`10^{22}`$ W cm-3.

The target explodes as a result of dynamic compression. At explosion an ultrabright quasipoint X–ray source is formed with registered energy range from $`2\mathrm{}3`$ keV up to 10 MeV and more, with maximum intensity in the range $`10\mathrm{}30`$ keV. Products of a target explosion precipitate on the internal walls of concentrator and on the special accumulating screen — metallic disk 1 cm in diameter. The precipitated products have the form of irregularly scattered drops, beads, films, etc.

### I.2 Analysis

Precipitated products were analyzed to determine their isotopic and elemental compositions. The analytical methods employed were:

* electron probe microanalysis (EPMA);
* Auger–electron spectroscopy (AES);
* secondary ion mass–spectrometry (SIMS);
* laser mass–spectrometry (LMS);
* Rutherford backscattering of accelerated alpha–particles (RBS).

As a rule, after the experiment, elemental analyzes of the surface and surface layer up to 3 $`\mu `$m deep of each accumulating screen (made of chemically pure material: Cu (99.99 % pure), Ag (99.99 %), Ta (99.97 %), Pb (99.75 %), etc.) were performed. Minimal number of microprobes was from 10 to 20, some regions (spots) were also analyzed by several methods.

By March 12, 2003, altogether 14520 analyzes were performed and registered. The summary of the analytical methods and number of analyzes are listed in Table 1.

### I.3 Main results of elemental composition analyzes after one of the typical experiments

Analyzes show that the products of impact explosion precipitated on the accumulating screen contain a wide spectrum of the light, medium and heavy chemical elements that are absent in comparative quantities or concentrations in the initial materials involved in the process of nuclear transformation.

This fact is illustrated in Table 2, where the quantities of impurities in chemically pure copper (used for the target and accumulating screen) are compared with those detected in the products precipitated on the accumulating screen after one of the typical experiments. Correlation ratios between the concentration of impurities in the target material and their quantities in the precipitated products have near zero or negative values.

The composition of the initial environment and target materials were measured with a mass–spectrometer of glow discharge VG 9000 (Thermo Elemental), with sensitivity from 100 % to ppt level in a single analysis.

Three different instruments were used for analysis of the elemental composition of the accumulating screen surface layer after the target impact explosion:

* Auger spectroscope JAMP–10S (JEOL, Japan) — 21 spots $`<`$1.0 $`\mu `$m in diameter and 0.002 $`\mu `$m deep were analyzed;
* EPMA analyzer REMMA 102 (Ukraine) — 113 areas 11$`\times `$11 $`\mu `$m and 3.0 $`\mu `$m deep, and 417 extrinsic beads of different shape on the surface of the accumulating screen.
* SIMS analyzer IMS 4f (CAMECA, France) — 5 areas 250$`\times `$250 $`\mu `$m and 0.4 $`\mu `$m deep (because SIMS is the destructive method, it was used after Auger and EPMA);

As can be seen from Table 2, both the light elements (with the mass number less than the initial copper) and the heavy ones are present. It should be emphasized that the elements with atomic masses exceeding two masses of the initial element (elements with mass number $`Z>58`$) were detected as well. Quantities of the mentioned elements by a few orders of magnitude exceed impurities in the initial material.

To get quantitative estimations for the number and the concentration of deposited atoms, the statistical analysis was performed for the following layers:

1. Thin surface layer about 0.002 $`\mu `$m;
2. Intermediate layer from 0.002 $`\mu `$m to 0.4 $`\mu `$m;
3. Rest of the material — layer from 0.4 $`\mu `$m to 3.0 $`\mu `$m (maximum depth of analyzes).

As the result of statistical evaluation, the concentrations in the above–mentioned layers were obtained for each element. The validation criterion was the accordance between the calculated element concentration and concentration measured by all applied instruments.

Statistical analysis of the quantity of deposited atoms revealed the following:

* Deposited atoms ($`2.8\times 10^{16}`$ atoms) concentration is maximum in the surface layer 0.002 $`\mu `$m thick (1st layer) and the concentration of copper (main target and screen material) is lower than 10 %.
* Main quantity of deposited elements ($`1.3\times 10^{18}`$ atoms) and the majority of the elements with masses exceeding two masses of the initial material (i.e. Cu) are present in the 2nd layer (From 0.002 to 0.25 $`\mu `$m).
* Quantity of deposited atoms in the 3rd layer (0.25 $`\mu `$m and down to the limit of the analyzer, 3.0 $`\mu `$m) is $`3.1\times 10^{17}`$ atoms.

Total number of deposited atoms, excluding copper, present in the substance up to 3.0 $`\mu `$m deep was calculated as sum of the atom quantities in the volumes of each layer and equals $`1.6\times 10^{18}`$ atoms. Change of the thickness of the 2nd layer containing the main part of deposited elements from 0.2 $`\mu `$m to 0.4 $`\mu `$m gives the range $`(1.4\mathrm{}2.4)\times 10^{18}`$ for the total number of deposited atoms.
* Quantity of deposited atoms in randomly dispersed particles on the surface of the accumulating screen is equal to $`6.0\times 10^{17}`$ atoms.

The total number of deposited atoms, excluding copper, after the experiment in comparison to the initial concentration is shown in Table 2.

The same estimation for the synthesized atoms was derived from the ‘marked target’ experiments. The targets for these experiments were fabricated from radioactive cobalt (60Co). After an impact, the system activity (number of radioactive decays per second) decreased by the amount equal to the transmutation of $`10^{18}`$ atoms of the active target volume — intensity of the 60Co spectral lines decreased while no other radioactive isotopes line did appear.

### I.4 Isotopic composition analyzes
of accumulating screens

To confirm that we deal with the nucleosynthesis process, the isotopic composition of the products on the surface of the accumulating screens was analyzed. The stimulus for this was the well–known fact that the isotopic composition of any element is almost identical throughout the Solar system (see, e.g. iupac ; wasser ).

Analyzes of the isotopic composition were done with two different methods, i.e., LMS and SIMS (see Table 1). As was found out, most of the analyzed spots of an accumulating screen had an isotopic composition that differs from a natural one. Examples of the isotopic composition for some elements are presented in Fig. 1.

As can be seen from this figure, the isotopic composition of the elements synthesized from the copper target, differs significantly from the natural one. Hence, this proves an artificial origin of the detected nuclei.

## II Synthesis of the elements at the border of the periodic table and beyond it

The presence of heavy elements (tantalum, tungsten, gold, and lead) among the products of nucleosynthesis from the copper target encouraged us to carry out the experiments with heavier targets. In fact, since lead ($`A(\mathrm{Pb})207`$) was produced from the copper target ($`A(\mathrm{Cu})64`$), one should expect in similar way to get elements at the border of the periodic table and even beyond it if a target made from heavier elements is used.

For this aim, the experiments with the platinum, lead, and bismuth targets were carried out. The results of analyzes revealed the presence of long-living isotopes of chemical elements at the border and even beyond the known part of periodic table.

To identify the elemental composition of the accumulating screen surface, the Auger spectroscopy (AES) method was used. The method is non-destructive and highly localized ($`50\mathrm{}100`$ nm). The depth of analyzes is low ($`1\mathrm{}2`$ nm). The method ensures a wide range of detectable elements (all but H and He) and a relatively high sensitivity ($`0.1\mathrm{}1`$ at.%).

The spectra obtained with the help of a Auger–microprobe JAMP–10S (JEOL, Japan) typically contained Auger peaks of a considerable amount of chemical elements that hindered their identification. To identify the Auger peaks, the measurement was performed in a wide range of energies ($`30\mathrm{}3000`$ eV, energy resolution $`0.5\mathrm{}1.2`$ %) to cover the maximum number of Auger peaks’ series, and prolonged exposures (up to 3 hours) were used to identify the low intensity peaks. The spectra were recorded in differential form at the accelerating voltage of the electron probe of 10 kV and beam current $`10^6\mathrm{}10^8`$ A. Residual pressure in the sample chamber was $`5\times 10^7`$ Pa. Artifacts like electrical charging, characteristic energy losses, and chemical shift were taken into consideration. A standard concentration calculation program supplied by the manufacturer (JEOL) was used for quantitative analysis.

In the process of analysis, the unidentified peaks with energies of 172, 527, 1096 eV, and a doublet of 130 and 115 eV were registered on the surface of accumulating screens. These peaks do not correspond to any of the catalogued peaks of chemical elements and cannot be referred to any of the known artifacts. One of the possible explanation of the mentioned peaks is the presence of long-living transuranium elements.

Two sections of the Auger spectrum (energies 172, 527 eV) out of six unidentified peaks are presented in Fig. 2.

To provide elemental analysis of the accumulating screens surface in a wide range of masses, a SIMS analyzer IMS 4f (CAMECA, France) was used.

Here are some key characteristics of the instrument:

* mass detection range — up to 480;
* mass resolution ($`M/\mathrm{\Delta }M`$) — 2500;
* locality — 0.2 $`\mu `$m;
* primary ions: Cs+, O$`{}_{2}{}^{}{}_{}{}^{+}`$, O-, Ar+;
* secondary beam stability (20 min): $`\mathrm{\Delta }I/I=0.7\%`$,
$`\mathrm{\Delta }M/M=7\times 10^6`$.

A complex approach was used for elemental analyzes of the accumulating screens surface and included studies of the topography of chemical elements on the surface in the cluster suppression mode (Offset) and high mass resolution mode simultaneously with the photo film recording of the mass distribution.

In the Offset mode, an additional voltage is applied and, owing to the potential barrier, allows us to identify monoatomic ions even in the low–mass resolution mode and, hence, to exclude clusters and multicharged ions.

From the whole mass–spectrum, we consider as unidentified only those peaks which (a) were not observed in the reference material spectrum, (b) did not coincide with the other elements’ distribution on the screen surface, (c) offset voltage increase led to an increase of the intensity in comparison with clusters’ peaks.

Typical mass spectrum in the mass range $`250\mathrm{}280`$ obtained by IMS 4f is presented in Fig. 3. Filled peaks, corresponding to 265 and 267, are unidentified; 268, 270, 272, 274 are Cu4O. Increasing the offset voltage led to a significant decrease of the intensity for copper clusters and the negligible fall of the intensity of 265 and 267 peaks indicating that they belong to monoatomic ions. To exclude possible chemical elements combinations, the photos of all present elements distributions were done, and they show that none of the unidentified masses coincide with any of them.

Up to date, 20 samples were thoroughly checked for unidentified masses. On 17 samples, more than 100 unidentified masses were found in the range $`221\mathrm{}475`$. Most frequently were detected the peaks of 271, 272, 330, 341, 343 masses. Accumulating screens that were used repeatedly in experiments with copper targets also revealed the presence of long-living isotopes in the mass range $`300\mathrm{}481`$.

In addition to above–mentioned spectrometric methods, the samples were examined for the presence of superheavy elements with the most direct method — Rutherford backscattering. The samples were irradiated with the beam of 27.2 MeV alpha–particles accelerated in a U–120 cyclotron Shvedov . The energy spectrum of scattered projectiles exhibited the scattering centers corresponding to mass numbers in the range $`200\mathrm{}400`$ (Fig. 4).

Thus, different methods of analysis detected the presence of unidentified superheavy elements among the products of artificial nucleosynthesis.

## III Summary of results and discussion

We have presented the survey of experimental studies of nucleosynthesis products obtained in experiments with the supercompressed substance. The essence of results is as follows:

* In substance compressed up to a superhigh density (presumably, $`10^{26}`$ cm-3), the process of nucleosynthesis and transmutation occurs. This process takes place over a macrovolume of the target substance ($`2\mathrm{}5\times 10^{18}`$ atoms/kJ).
* Nuclear processes in a target are collective and multiparticle that is proved by significant amount ($`>10^{13}`$) of synthesized nuclei with masses more than two masses of the initial nucleus.
* No $`\alpha `$–, $`\beta `$–, $`\gamma `$–active isotopes were observed in the products of laboratory nucleosynthesis, the radiation intensity never exceeded the background intensity.
* The activity of the targets marked with radioactive isotopes was reduced after compression impact by the value equal to the transmutation of $`10^{18}`$ atoms of the active target zone (collapse zone) for every 1 kJ of the driver energy.
* Developed installation has high reproducibility in reaching the conditions in a compressed substance necessary for the ignition of collective multiparticle fusion–fission reactions of the full spectrum of chemical elements.
* The products of nuclear transmutation reveal the presence of long–living isotopes of superheavy elements.

We consider that these reactions are, in some sense, similar to the pycnonuclear reactions that take place under the specific conditions of compact astrophysical objects (white dwarfs, supernovas, outer shells of neutron stars, etc.).

The new challenge is the development of a method and a technique for simultaneous measurement of both the atomic and nuclear properties of superheavy atoms, i.e., nucleus charge and mass.

The results of our experiments allow us to suppose that the synthesized nuclei are the product of clusterization in the decaying superdense electron–nucleon plasma that corresponds to maximization of the binding energy per nucleon, dependent on the substance density, on the one hand, and on the neutron concentration in the certain volume of clusterization, on the other hand.

## Acknowledgments

It is a pleasure for the authors to thank
E.V. Bulyak, I.A. Kossko, Y.V. Syten’ko,
S.A. Shestakov, D.G. Birukov, V.V. Kovylyaev,
M.I. Gorodyskiy, S.S. Ponomarev, G.A. Zykov


text:
# Efficient frequency doubler for the soft X-ray SASE FEL at the TESLA Test Facility

## 1 Introduction

The soft X-ray FEL at the TESLA Test Facility at DESY (TTF FEL) will cover a spectral range between approximately 60 nm and 6 nm wavelength. The minimum wavelength of 6 nm is determined by the maximum electron beam energy of 1 GeV. It would be extremely interesting to extend this range into the so-called water window, i.e. the range between the K-Absorption edges of carbon and oxygen at 4.38nm and 2.34 nm, respectively. This would allow time-resolved studies of organic molecules above the carbon K-edge, and, for example, the investigation of thick, hydrated biological samples by X-ray microscopy without the need of shock-freezing . There are several techniques under consideration to reach this goal: A X-ray SASE FEL at shorter wavelengths, which requires higher beam energy or a shorter undulator period, and the generation of harmonics through a nonlinear mechanism driven by bunching at the fundamental.

There are basically two methods to up-convert the fundamental radiation frequency via nonlinear harmonics. These methods are the generation of third harmonic radiation in a planar SASE undulator through a nonlinear mechanism driven by bunching at the fundamental frequency, and a two-undulator (second) harmonic generation scheme, also referred to as the ”after-burner” method.

SASE FELs are capable to produce powerful radiation not only at the fundamental frequency, but also at higher harmonics. When a beam is strongly bunched in the sinusoidal ponderomotive potential formed by the undulator field and the radiation field of the fundamental frequency, the electron beam density spectrum develops rich harmonic contents. Coherent radiation at the odd harmonics can be generated in a planar undulator and significant power levels for the third harmonic can be reached before the FEL saturates (see  and references therein). It is expected that the power of the transversely coherent third-harmonic radiation can approach the 1% level of the fundamental power at the TTF. Since the nonlinear harmonic generation occurs naturally in a planar undulator, no additional FEL hardware components are required for this method.

An idea of using two undulators, with the second undulator resonant to one of the harmonics of the first one, was considered in  (it is also referred to as the “after-burner” method). The first undulator is long enough to reach saturation and produce strong spatial bunching in harmonics. The bunched beam generates coherent radiation in the second undulator which follows immediately the first one. The main problem with this approach is the large induced energy spread which will be generated by the bunching of the electron beam at the fundamental frequency. While this energy spread is necessary for the bunching, it degrades the performance of the radiator section at the harmonic frequency. Another method to generate higher harmonics is the high-gain harmonic generation scheme (see  and references therein). It should be noted that all previous studies of such frequency multiplication schemes were performed in the framework of the steady-state approach, and the question arises how effectively they work in the case of a SASE FEL.

In this paper we propose an effective frequency doubler scheme for a SASE FEL. It consists of an undulator tuned to the first harmonic, a dispersion section, and a tapered undulator tuned to the second harmonic (see Fig. 1). The first stage is a conventional soft X-ray SASE FEL. The gain of the first stage is controlled in such a way that the maximum energy modulation of the electron beam at the FEL exit is about equal to the local energy spread, but still far away from saturation. When the electron bunch passes through the dispersion section this energy modulation leads to effective compression of the particles. Then the bunched electron beam enters a tapered undulator, and from the very beginning produces strong radiation because of the large spatial bunching. The strong radiation field produces a ponderomotive well which is deep enough to trap the particles, since the original beam is relatively cold. The radiation produced by these captured particles increases the depth of the ponderomotive well, and they are effectively decelerated. As a result, much higher power can be achieved than for the case of a uniform undulator. In addition, the output radiation exhibits excellent spectral properties due to the suppression of the sideband growth in a tapered undulator.

In this paper we analyze for the first time the frequency multiplication scheme for a SASE FEL. Simulations using the time-dependent code FAST , upgraded for the simulation of higher harmonics, provide a “full physics” description of the process. We illustrate the operation of the proposed frequency doubler for the parameters of the TTF FEL. The result of our study is that a frequency doubler can be implemented at the TESLA Test Facility. With an additional, 13.5 m long undulator (with 1.95 cm period and 0.39 T peak magnetic field), the TTF FEL would be able to produce radiation down to 3 nm wavelength with an output peak power in the GW range, and with excellent spectral properties. The power of the third harmonic from this device (i.e. at 1 nm wavelength) is still in the ten MW range, exceeding any other pulsed radiation source presently available at 1 nm, and sufficiently high for novel applications. In fact, this device approaches the operating range of the future TESLA XFEL.

The paper is organized as follows. Section 2 summarizes recent updates of the TTF FEL parameters which refer mainly to a smaller local energy spread than it was expected before. Section 3 describes the operation of the frequency doubler scheme. Section 4 discusses the integration of the frequency doubler scheme into the TESLA Test Facility. Appendixes 1 and 2 contain information on alternative ways for attaining 3 nm wavelength at the TTF FEL using a uniform undulator and an after-burner, respectively.

## 2 6 nm operation of the SASE FEL at DESY

A soft X-ray SASE FEL will be commissioned as a user facility in 2004 at the TESLA Test Facility at DESY, covering the wavelength range down to 6 nm. The first description of the TESLA Test Facility FEL has been written in 1995 . An update of the TTF FEL design was fixed in 2002 . The main issues were the change of the undulator focusing structure (separated instead of integrated) and the bunch compression system (removing the low-energy bunch compressor and introducing a RF structure operating at the 3rd harmonic).

The analysis of recent experimental results obtained at the TTF FEL, Phase 1  led to the idea that the value of the local energy spread is significantly less than it was expected originally. Only a small value of the local energy spread allows one to obtain a high-peak value of the bunch current after a single bunch compressor. Direct measurements of the local energy spread gave an upper estimate below 5 keV at a beam current of about 100 A . An extrapolation of this experimental result to the project value of 2500 A beam current gives an estimate for the local energy spread of about 100 keV. This is in agreement with recent start-to-end simulations of the bunch in the TTF linac  predicting the value of the local energy spread within the lasing part of the bunch to be well below 200 keV (see Fig. 3). Note that the original project value for the local energy spread was 1 MeV . The local energy spread is one of the essential parameters influencing the FEL operation, and we have to analyze which impact it may have on the operation of the current version of the TTF FEL and its future developments.

Tables 1 and 2 show an updated list of parameters for the SASE FEL at the TESLA Test Facility at DESY. The decrease of the local energy spread results in shorter a saturation length as it is seen in Fig. 5. This is a consequence of the fact that the previous design did not have sufficient safety margin with respect to the energy spread. Figure 5 shows the time and spectral structure of the radiation pulse for the linear and saturation regime.

With a small value of the local energy spread we fall in a different region of physical parameters which reveals a possibility to implement different FEL amplifier schemes. In particular, the application of dispersion sections for beam bunching becomes very effective.

## 3 Operation of a frequency doubler

In this section we illustrate the operation of the frequency doubler for the parameters of the TTF FEL, Phase 2. The beam parameters are given in Table 1, and the parameters of the doubler undulator and the dispersion section are listed in Table 3. The resonance wavelength in the main X-ray undulator is equal to 6 nm, and the doubler undulator is tuned to the resonance wavelength of 3 nm.

The frequency doubler scheme operates as follows. The electron bunch enters the main X-ray undulator and produces SASE radiation at the wavelength of 6 nm. During the amplification process the radiation power grows exponentially with the undulator length. Simultaneously, the energy and density modulation of the beam are growing. At the end of the X-ray undulator the beam energy modulation is comparable with the local energy spread. In the present example this takes place at an undulator length of 12.3 m. The output characteristics of the radiation are illustrated with plots in Fig. 5 (upper row). The upper left plot in Fig. 6 shows the phase space distribution of particles for the slice corresponding to 330 fs time coordinate along the bunch. Such a picture is typical for every spike. We see that the modulation amplitude at the second harmonic is small (see lower left plot in this figure). On the other hand, there is visible energy modulation with an amplitude of about the value of the local energy spread. When the electron bunch passes through the dispersion section this energy modulation leads to effective compression of the particles as it is illustrated with plots in the right column in Fig. 6. When the bunched beam enters the undulator tuned to the second harmonic, it immediately starts to produce powerful radiation at the second harmonic. That is important, but not the main feature of our proposal. Let us study more closely the phase space distribution of the particles at the exit of the dispersion section. We see that this distribution is double-periodic with respect to the second harmonic. I.e., only each second bucket is populated with particles. If we trace the evolution of the FEL process in the uniform undulator, we find that “thermalization” takes place. The population of particles in the originally filled buckets is reduced, and some of these particles travel into the ponderomotive well of originally empty buckets. Simultaneously, the energy spread in the bunch grows due to the FEL process. Finally, saturation occurs, and the properties of the radiation are close to those described in Appendix 1.

Thus, another key idea of our proposal is to preserve the original bunching of the beam and to organize an effective extraction of the energy from the electron beam. This idea is simply realized by using a tapered undulator. The process of amplification proceeds as follows. The bunched beam (double-periodic) produces strong radiation from the very beginning because of the large spatial bunching. The strong radiation field produces a ponderomotive well which is deep enough to trap particles, since the original beam is relatively cold. The radiation produced by these captured particles increases the depth of the ponderomotive well, and they are effectively decelerated. The undulator tapering preserves the synchronism of trapped particles and radiation, and a significant fraction of the energy can be extracted. Simulations using the time-dependent FEL code FAST  confirm this qualitative picture. Figure 8 shows the evolution of the energy in the radiation pulse in the tapered undulator. Despite of the original spiking seeding the process of the second harmonic (see Fig. 6), we effectively trap a significant fraction of the particles, and can achieve much higher power than for the case of an untapered undulator (see Appendix 1). Figure 8 shows the temporal and spectral structure of the radiation pulse after 9 and 13.5 meter of second-harmonic undulator. One can see that the radiation pulse length is about a factor of two shorter than that of a traditional SASE FEL (see Fig. 10). This is a consequence of the nonlinear transformation provided by the dispersion section. Another important feature of the radiation from a tapered undulator is the significant suppression of the sideband growth in the nonlinear regime (compare Figs. 8 and 10). This means that in the proposed scheme the spectral brightness of the radiation is increased proportionally to the radiation power. In the case of a uniform undulator the peak brightness is reached at the saturation point and is then reduced due to the sideband growth .

## 4 Implementation of a frequency doubler scheme at the TESLA TEST Facility at DESY

The TTF FEL is a pioneer user facility in the VUV and soft X-ray wavelength range. The conceptual design was first elaborated in 1995 . At that time only a single-pass SASE FEL was considered. In 1998 it was decided to upgrade the TTF FEL by a seeding option . Other possible extensions of the TTF FEL user facility making use of some space left behind the main FEL undulator have been discussed .

This free space downstream of the main undulator (see Fig 2) would be sufficient for a 9 meter long ($`2\times 4.5`$ m sections) second harmonic radiator. It is planned now that the main X-ray undulator will consist of six sections. This is mainly dictated by the wish to have sufficient safety margin during the first commissioning of the SASE FEL. However, our analysis in the present paper shows that for the expected electron beam parameters five undulator modules will be more than sufficient. This means that in the future there will be space for three modules of the frequency doubler. This will allow to realize the full frequency doubler and to extend the wavelength range of the TTF FEL down to 3 nm with GW-level power without changing the present TTF layout. Additional hardware to be manufactured is a dispersion section and a 13.5 m long undulator. Since the required strength of the dispersion section is small, its design may be similar to that of a phase shifter for the X-ray SASE FEL . It will also be necessary to include a linear tapering of 0.14%/m ( in the 13.5 m long undulator with 1.95 cm period and 0.39 T peak magnetic field. This does not appear to be a serious problem. In conclusion we can state that these are minor expenses in view of the significant extension of the capabilities of the TTF FEL user facility.

## 5 Conclusion

In this paper we have described an effective frequency doubler for SASE FELs. For the first time the frequency multiplication scheme has been analyzed for SASE FELs. To be specific, we have illustrated the proposed FEL scheme using the parameters of the TESLA Test Facility. It is shown that a frequency doubler allows to reach the water window and to produce GW-level radiation pulses. It is important to note that a frequency doubler would be well compatible with the seeding option, and would transfer all the advantages of the seeding option to twice the photon energy . In general, the frequency doubler scheme has several significant advantages over traditional SASE FELs with uniform undulators:

* shorter total magnetic length;
* a possibility to attain higher output power and brightness of the radiation;
* shorter radiation pulse duration.

The realization of a frequency doubler at the TTF is considered as extremely important not only for reaching the water window. It may also be relevant for realizing the X-ray FEL user facility in the 0.1 nm range. The use of frequency doublers at a large-scale facility could result in a significant reduction of the project costs. In addition, the safety margin for facility operation becomes more relaxed.

## Acknowledgments

We thank J. Krzywinski, J. Rossbach, and M. Tischer for many useful discussions. We thank C. Pagani, J.R. Schneider, and D. Trines for their interest in this work.

## Appendix 1: 3 nm option with uniform undulator

The maximum energy of the TTF linac is limited to 1 GeV. The present TTF undulator has a period of 2.73 cm and a peak magnetic field of 0.47 T which limits the minimum wavelength to 6 nm. In principle, one can consider a scenario of constructing a short-period undulator allowing to generate wavelengths down to 3 nm. Such an option was not studied previously because of the large value of 1 MeV for the expected local energy spread. However, in view of the change of this parameter it seems to be reasonable to perform such a study in order to have a possibility for comparison of different approaches for attaining 3 nm wavelength. This section illustrates the option of a 3 nm SASE FEL with a uniform undulator. The parameters of the electron beam are presented in Table 1, and those of the undulator are given in Table 4.

Figure 10 shows the evolution of the energy in the radiation pulse along the undulator. The growth of the radiation energy starts to saturate at an undulator length of about 36 m, and the radiation energy is about 100 $`\mu `$J. Further increase of the undulator length leads to a slow increase of the radiation energy, while the spectral width is increased due to sideband growth in the nonlinear regime (see Fig. 10).

## Appendix 2: 3 nm option with after-burner

An option of second harmonic generation at TTF FEL using an after-burner has been under discussion for several years . Here we present the results of time-dependent simulations of the after-burner scheme. The parameters of the electron beam are given in Table 1, and the parameters of the undulator and the output radiation are presented in Table 5.

In the after-burner scheme the spent electron beam leaving the main X-ray undulator passes an undulator tuned to the second harmonic. At the exit of the main undulator the electron beam has a pronounced amplitude of density modulation at the second harmonic which serves as input signal for the second-harmonic undulator (see Fig. 13). When the electron beam enters the after-burner radiator, it readily starts to produce radiation. However, the power growth saturates quickly at 2 meters, as it is seen from Fig. 13. This is due to the large energy spread induced in the main undulator (see Fig. 13). Figure 13 presents the temporal and spectral structure of the radiation pulse at the exit of the after-burner. We find that the FWHM pulse duration is relatively large, about 250 fs. The spectral width is also large, about 0.4%, and is driven by the FEL process in the main undulator. Finally, the level of output radiation power is low, about 60 MW, i.e. about 3% of the power of the fundamental frequency of the main undulator. Remembering that SASE radiation from a planar undulator at saturation always contains the 3rd harmonic at a per cent level, we conclude that harmonic generation using an after-burner does not provide extra opportunities for user applications.


text:
# Dynamical Properties of the Delta Kicked Harmonic Oscillator

## Abstract

We propose an efficient procedure for numerically evolving the quantum dynamics of Delta Kicked Harmonic Oscillator. The method allows for longer and more accurate simulations of the system as well as a simple procedure for calculating the system’s Floquet eigenstates and quasi-energies. The method is used to examine the dynamical behaviour of the system in cases where the ratio of the kicking frequency to the system’s natural frequency is both rational and irrational.

Kicked Hamiltonian systems are firmly established as prototypical models for studying chaos. Enormous effort has been spent examining models like the standard map and the kicked rotator. The study of kicked degenerate systems, those that do not obey the KAM theorem, have received relatively little attention but display fascinating dynamical properties both classically and quantum mechanically. The Kicked Harmonic Oscillator (KHO) is an example of such a system. The system was originally proposed as a 2-dimensional model of kicked charges in a uniform magnetic field Zaslavsky . It has also been proposed as a model for electronic transport in semiconductor superlattices Fromhold  and atom optic modeling using ion-traps Gardiner . In particular, the quantum kicked harmonic oscillator, in the same way as the quantum kicked rotator Levi , may be simulated efficiently on a quantum information processor using the Quantum Fractional Fourier Transform Klappenecker .

The Hamiltonian of the system may be written as

$$H_{KHO}=\omega _{HO}\left(\frac{p^2}{2}+\frac{q^2}{2}\right)+V(q)\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta (tnT_K)$$
(1)
where the kicking potential $`V(q)`$ in our system is given by

$$V(q)=\mu \mathrm{cos}(kq)$$
(2)
where k is a parameter related to $`\omega _{HO}`$ and the mass $`m`$. It is distinct from other widely studied systems because of the presence of two different frequencies in the Hamiltonian, namely the natural frequency of the Harmonic Oscillator $`\omega _{HO}`$, and the kicking frequency $`\omega _K=2\pi /T_K`$. It is the ratio between these frequencies, $`\omega _{HO}/\omega _K=1/R`$, and the strength of the kicking potential, that determines the global behavior of the system Zaslavsky . It is appropriate to note here that the kicking potential we use is even. Using an odd kicking potential, say $`V(q)=\mu \mathrm{sin}(kq)`$, changes the global behaviour of the system dramatically Dana2 . We do not attempt to analyse this case although the numerical techniques that we present later can be applied to any delta kicking potential.

Understanding the quantum system analytically has proved to be problematic. While expressions have been found for the Floquet operator Daly , these generally involve complicated infinite summations of Bessel functions and are too cumbersome to yield any important information about the dynamics of the system. Numerical simulations are therefore essential to understanding how the system evolves. However, in general, factors such as memory and computation time reduce the accuracy of simulations before much can be concluded about the eigensolutions or long time dynamical evolution of the system. We outline a procedure, based on the Fractional Fourier Transform, that allows for efficient numerical simulation of the system’s evolution as well as an effective and straightforward approach to calculating the system’s Floquet operator. We use this approach to confirm previous results for the system where $`R=4,5`$ as well as showing that delocalisation can occur when $`R`$ is irrational if the kicking potential is strong enough.

The Floquet operator for this system can be factorized into two parts:

$$U_{KHO}(T_K)=U_{SHO}(T_K)U_{kick}$$
(3)
where because of the delta function in Eq.(1) we can express the kick operator as

$$U_{kick}=\mathrm{exp}\left(\frac{i}{\mathrm{}}\mu \mathrm{cos}(kq)\right),$$
(4)
which is clearly diagonal in the position basis and between kicks

$$U_{SHO}(T_K)=\mathrm{exp}\left(\frac{i\omega _{HO}T_K(p^2+q^2)}{2\mathrm{}}\right),$$
(5)
which is diagonal in the Fock state basis but not in the position basis or the momentum basis. Changing from the number basis to the position basis requires a matrix-vector multiplication ($`N^2`$ operations) where $`N`$ is the number of dimensions being used to approximate the infinite dimensional Hilbert space. However, it is well known that changing between the position and momentum basis can be done through Fourier Transforms. This fact is used to great advantage in the split step method Feit . In this method the exact operator $`U_{SHO}`$,applied for the small duration $`\mathrm{\Delta }t`$, can be approximated by

$$U_{SHO}^{}(\mathrm{\Delta }t)=\mathrm{exp}\left(\frac{i\omega _{HO}\mathrm{\Delta }tp^2}{4\mathrm{}}\right)\mathrm{exp}\left(\frac{i\omega _{HO}\mathrm{\Delta }tq^2}{2\mathrm{}}\right)\mathrm{exp}\left(\frac{i\omega _{HO}\mathrm{\Delta }tp^2}{4\mathrm{}}\right)$$
(6)
where we use the Fast Fourier Transform (FFT) Cooley , to perform the change of basis. Although this method is not exact, the error is $`O(\mathrm{\Delta }t)^3`$, and can therefore be made acceptably small by an appropriate choice of $`\mathrm{\Delta }t`$. The number of operations needed to evolve the system from time $`t`$ to time $`t+\mathrm{\Delta }t`$ is $`O(N\mathrm{log}N)`$. However, using the stability condition for the Split-Step method Weideman , $`\mathrm{\Delta }t<\frac{l^2}{\pi }`$, where $`l`$ is the grid spacing, we roughly estimate the number of operations needed to accurately evolve the system over a finite time to be of $`O(N^2\mathrm{log}N)`$. The essential point is that we avoid the use of matrix-vector multiplications through the use of Fast Fourier Transforms. This technique is relatively efficient computationally and does not require large amounts of storage.

It is possible however, to exploit certain symmetries in the Hamiltonian to evolve the system with far greater efficiency. The previous discussion regarded the Fourier Transform as a passive transformation, that is, one that changes the position basis to the momentum basis. However, for the harmonic oscillator, a Fourier Transform can also be regarded as an active transformation. In 1980, Namias defined the Fractional Fourier Transform in terms of the Hermite polynomials and showed that it is in fact the time evolution operator of the harmonic oscillator up to a phase Namias . Explicitly, taking $`\mathrm{\Psi }(x,t)`$ as the time-dependent wavefunction of the harmonic oscilator, he showed that

$$\mathrm{\Psi }(x,t)=U_{SHO}(t)\mathrm{\Psi }(x,0)=e^{it/2}F_t\mathrm{\Psi }(x,0),$$
(7)
where the operator $`F_t`$ is a generalisation of the ordinary Fourier transform Defn . Indeed the latter is just the special case when $`t=\pi /2`$. This result is important for our system because of the existence of Fast Fractional Fourier Transforms (FFrFT’s) Ozaktas . There are a number of algorithms in the literature for performing the Fast Fractional Fourier Transform as efficiently as the regular FFT code . This means that, in the position basis, the evolution of the wavefunction over a finite time may be written as

$$\mathrm{\Psi }(x,T_K)U_{KHO}(T_K)\mathrm{\Psi }(x,0)=e^{iV(x)/\mathrm{}}e^{iT_K/2}F_{T_K}\mathrm{\Psi }(x,0),$$
(8)
and can be performed in about $`O(N\mathrm{log}N)`$ operations. This allows us to perform highly efficient numerical simulations of pure state quantum dynamics in Hilbert space dimensions of about $`2^{17}`$. This enables us to follow the diffusing quantum evolution for larger times and in greater detail. Simulations of such a size would be prohibitively costly in memory and processor time using other simulation techniques.

As an additional advantage, the Floquet matrix itself can be calculated. This is a consequence of Eq. (8) and the fact that the Fractional Fourier Transform matrix can be generated without much additional overhead Ozaktas . This allows for quick and straightforward calculation of the Floquet eigenstates and eigenvalues using the standard LAPACK routines.

One of the interesting properties of the system given in (1) and (2), is the presence of a stochastic web that spans all of the classical phase space for certain values of the parameter $`R`$. Indeed the values of $`R`$ for which this crystalline structure appears can be seen to be related to the tessellation of the phase plane, where the values of $`R=3,4,6`$ represent filling (tiling) the plane with triangles, squares or hexagons respectively (we ignore the degenerate cases of $`R=1,2`$). In the quantum mechanical problem these specific values of $`R`$ have, for all values of the kicking strength, been shown to give rise to extended quantum eigenstates and therefore delocalisation of the quantum state over time Berman ; Borgonovi ; Frasca .

In the other cases where $`R`$ is an integer, a quasi crystalline structure is observed in the classical phase plane.In the quantum system there is a large suppression of the states quantum evolution with respect to energy Borgonovi ; Shepelyansky . These results do not appear to contradict those of the tight binding prediction Frasca , where extended states are predicted for all rational values of the frequency ratio. Indeed, there is a value of the kick strength, above which, the quantum system diffuses at a rate, less than, but comparable to that of the classical system Shepelyansky , whose energy is known to grow linearly with time.

In the case of irrational frequency ratios all traces of a crystalline structure disappear in the classical phase plane. The tight binding model predicts localisation for all values of the kick strength Frasca , and this has been supported numerically for low kick strengths Borgonovi . As far as we are aware there is no published numerical data for irrational frequency ratios for this system with large kick strength.

To study the localisation properties of the system we first investigate the mean energy after each application of the Floquet operator. We calculate

$$E(t)=\psi _t|\frac{p^2}{2}+\frac{q^2}{2}|\psi _t,$$
(9)
as a function of the discrete time $`T_K`$ or the kick number $`n`$. A direct comparison may be made between the quantum state and a classical ensemble. Defining $`\alpha =(q+ip)/\sqrt{2\mathrm{}}`$ we have the wavefunction of a coherent state centered at the coordinates $`(q,p)`$,

$$\psi (q)=\left(\frac{\omega _0}{\pi \mathrm{}}\right)^{1/4}\mathrm{exp}\left(\frac{\omega _0}{2\mathrm{}}q^2+\sqrt{\frac{2\omega _0}{\mathrm{}}}\alpha q\frac{|\alpha |^2}{2}\frac{\alpha ^2}{2}\right).$$
(10)
The Wigner function of this state is a 2-D Gaussian that is centered about the point $`(q,p)`$ Hillery . The Husimi function is very similar to this and can be regarded as the convolution of the Wigner function with a 2-D Gaussian centered at the origin. In the numerical simulations we compare the quantum evolution of a coherent state with the classical evolution of an ensemble of normally distributed particles initially centered around $`(q,p)`$ and possessing the same marginals as the corresponding Husimi function.

When $`R=4`$, the quantum system can be shown analytically to be identical to the well studied Kicked Harper Model (KHM) Dana2 . In this case the quantum system diffuses for all values of the kick strength $`\mu `$, and for all initial conditions. To demonstrate this we center the initial state on a classical stable point $`(q,p)=(0,0)`$, and note that the energy of the quantum state grows with time (or kick number $`n`$), while the energy of the classical state remains bounded. From Fig. 1(a), we see that the quantum system is able to tunnel through the separatrix which completely isolates classical trajectories. For completeness we also illustrate how both the quantum and classical systems diffuse when placed on an unstable point $`(q,p)=(0,\pi )`$, on the classical system’s stochastic web. It is worth noting that at $`\mu =0.5`$, the stochastic web is too small to allow classical diffusion. However we see that the quantum system uses the separatrix to significantly increase its rate of diffusion.

In the case where $`R=5`$, the classical phase plane shows a quasi-crystalline structure. From Fig. (2) we can see that at low kicking strengths diffusion of the quantum system is somewhat suppressed. However, above a certain critical value of the kick strength $`\mu `$, which depends on the initial position of the quantum state, this suppression no longer seems to exist.

We see a similar type of global behavior when the frequency ratio is irrational even though the structure of the phase space patterns are quite different. Both the quantum and classical systems remain bounded up to different critical values of the kick strength $`\mu `$ (again the precise value of this parameter depends on the initial conditions). Above these values the the system begins to diffuse linearly. In the same way as before, there exists a range of kick strengths where the classical system diffuses and the quantum motion remain suppressed,see Fig.3(b)(c). It would appear that, in this region, the quantum state cannot flow along the same chaotic trajectories as the classical ensemble. This is opposite to the case for which $`R=4`$, where at low kick strengths, the quantum system appears to exploit the infinitesimal stochastic web to aid its dispersion, see Fig.1(a)(b).

As we have mentioned, using the definition of the Floquet operator given in Eq. (8), we can easily find a finite dimensional matrix that can then be diagonalised by brute force. This method has allowed us find eigenvalues and eigenvectors of matrices up to 3000x3000. To go larger requires more RAM. This can then be used to examine the structure of these eigenstates in detail as well as examining the system’s eigenvalue statistics. As an example we show that, in the case of irrational kick frequencies, the Floquet operator contains some eigenstates that are clearly localised in space and some that appear to be extended. The existence of these states highlights what we have already found in the diffusion simulations. At low kick strength (say $`\mu =1`$), Fig. 4, an initial state may have eigenstate components which are localised. This is not the case however if we examine the eigenstates for the system with a larger kick strength. When $`\mu =6`$ even the lower energy eigenstates appear to have some sort of extended nature Fig. 5. These figures do not show all of the individual eigenstate. We have zoomed in on a small section to emphasis the Husimi function’s structural similarity to classical phase map.

In conclusion, we have addressed the numerical problems associated with studying the quantum delta kicked harmonic oscillator. We have shown that the Fractional Fourier Transform may be used to infer crucial numerical information about the system. This includes improving the speed and accuracy with which we can evolve the system as well as giving an easy and effective way to analyse the system’s Floquet eigenstates and quasi-energies. This improved numerical technique will assist in future studies of this systems dynamics using either a quantum information processor or the various physical systems mentioned in the introduction.

Acknowledgements: This research was supported by funding from Enterprise Ireland under the research grant EI SC/2001/059.