text:
# Computer-Generated Holographic Optical Tweezer Arrays

## I Introduction

Since their invention in 1986 , optical tweezers have become increasingly valuable tools for research in the biological  and physical  sciences. Using a focused beam of light to trap and move matter, optical tweezers offer convenient, non-invasive access to processes at the mesoscopic scale. Most applications, however, have involved manipulating small numbers of particles or small volumes of soft materials because existing optical tweezer implementations can create just a few tweezers at once. Were they readily available, large arrays of optical tweezers could be used to organize microscopic particles into complex structures, to sort them intelligently, to study collective behavior in many-body systems, and to manipulate materials too delicate to trap with a single tweezer. We recently described  a method to create arrays of optical tweezers using computer-generated holographic beam splitters. This Article further explains how to design and fabricate the necessary holograms and how to integrate them into *holographic optical tweezer arrays* capable of trapping hundreds of particles simultaneously.

## II Optics of Optical Tweezers

An optical tweezer traps particles with forces generated by optical intensity gradients. Dielectric particles polarized by the light’s electric field are drawn up the gradients to the brightest point. Reflecting, absorbing and low-dielectric particles, by contrast, are driven by radiation pressure to the darkest point. Optically generated forces strong enough to form a three-dimensional trap can be obtained by bringing a laser beam with an appropriately shaped wavefront to a tight focus with a high numerical aperture lens. Microscope objective lenses offer an ideal combination of minimal aberration and large numerical aperture and often serve as the focusing element in practical implementations of optical tweezers  and variants such as the optical vortex .

The challenge in constructing an optical tweezer is to direct a laser beam into the objective lens’ back aperture so that the beam fills the aperture and so that its axis coincides with the optical axis in the aperture’s plane, at the point labelled B in Fig. 1. If the beam follows the optical axis, then it comes to a focus and forms a trap in the center of the lens’ focal plane. If, on the other hand, it enters the back aperture at an angle, the resulting trap is offset from the center of the focal plane, as indicated schematically in Fig. 1.

Directing the beam into the objective with a dichroic mirror allows other wavelengths to pass through unimpeded and can be useful for imaging the trapped particles, as in Fig. 1. The problem remains, however, of aiming the beam.

The telescope formed by lenses L1 and L2 in Fig. 1 addresses this problem by creating a conjugate point, B, to the back aperture’s center, B, at a convenient location. A beam of light passing through B also passes through B and forms an optical trap. In our implementation, L1 and L2 are high quality plano-convex lenses with 250 mm focal lengths. Such long focal lengths help to minimize aberrations, particularly longitudinal spherical aberration, which would be detrimental to trapping . More compact optical trains would require additional attention to minimizing wavefront distortions. References  and  offer more detailed discussions of this aspect of the optical design.

Multiple beams passing through B all pass through B and thus all form optical tweezers. A diffractive optical element (DOE) at B, as shown in Fig. 1, can split a single collimated laser beam into any desired distribution of beams, each emanating from B at a different angle, and thus each forming a separate trap . Figure 1 shows the computer-generated pattern for a binary phase hologram together with a photomicrograph of colloidal particles trapped in the resulting array of optical tweezers. The remainder of this Article addresses the theory and practice of creating holograms such as the example in Fig. 1 suitable for projecting arbitrary arrangements of optical tweezers.

## III Holographic Tweezer Arrays

### A Fourier Optics

A planar array of optical tweezers can be described by the intensity distribution, $`I^f(\stackrel{}{\rho })`$, of laser light in the focal plane of a microscope’s objective lens. This pattern is determined by the electric field of light incident at its input plane, as depicted in Fig. 2. Suppose that the input plane is illuminated by monochromatic light of wavelength $`\lambda `$. Its wavefront at the input plane, $`E^{in}(\stackrel{}{r})`$, contains both phase and amplitude information,

$$E^{in}(\stackrel{}{r})=A^{in}(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }^{in}(\stackrel{}{r})],$$
(1)
where the amplitude, $`A^{in}(\stackrel{}{r})`$, and phase, $`\mathrm{\Phi }^{in}(\stackrel{}{r})`$, are real-valued functions. The electric field in the focal plane has a similar form,

$$E^f(\stackrel{}{\rho })=A^f(\stackrel{}{\rho })\mathrm{exp}[i\mathrm{\Phi }^f(\stackrel{}{\rho })],$$
(2)
so that $`I^f(\stackrel{}{\rho })=|E^f(\stackrel{}{\rho })|^2=|A^f(\stackrel{}{\rho })|^2`$. These fields are related by the Fourier transform pair

$`E^f(\stackrel{}{\rho })`$ $`={\displaystyle \frac{k}{2\pi f}}e^{i\theta (\stackrel{}{\rho })}{\displaystyle d^2rE^{in}(\stackrel{}{r})e^{ik\stackrel{}{r}\stackrel{}{\rho }/f}}`$ (3)
$`\{E^{in}(\stackrel{}{r})\},\mathrm{and}`$ (4)
$`E^{in}(\stackrel{}{r})`$ $`={\displaystyle \frac{k}{2\pi f}}{\displaystyle d^2\rho e^{i\theta (\stackrel{}{\rho })}E^f(\stackrel{}{\rho })e^{ik\stackrel{}{r}\stackrel{}{\rho }/f}}`$ (5)
$`^1\{E^f(\stackrel{}{\rho })\},`$ (6)
where $`f`$ is the focal length of the lens and $`k=2\pi /\lambda `$ is the wavenumber of the incident light. The additional phase profile, $`\theta (\stackrel{}{\rho })`$, due to the lens’ geometry does not contribute to $`I^f(\stackrel{}{\rho })`$ and may be ignored without loss of generality .

### B Phase-Only Holograms

Obtaining a desired wavefront in the focal plane requires introducing the appropriate wavefront in the input plane. Most lasers, however, provide only a fixed wavefront,

$$E_0(\stackrel{}{r})=A_0(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }_0(\stackrel{}{r})].$$
(7)
Shaping $`E_0(\stackrel{}{r})`$ into $`E^{in}(\stackrel{}{r})`$ involves modifying both the amplitude and phase at the input plane. Changing the amplitude with a passive optical element necessarily diverts power from the beam and diminishes trapping efficiency. Fortunately, optical trapping relies on the beam’s intensity and not on its phase. We can exploit this redundancy by setting $`A^{in}(\stackrel{}{r})=A_0(\stackrel{}{r})`$ and modulating only the phase of the input beam to obtain the desired trapping configuration.

Several techniques are available for achieving the necessary phase modulation, and some of the associated practical considerations are discussed in Section VI. For the purposes of the present discussion, we will refer to the phase modulating element as a hologram or a diffractive optical element and treat it as if it acts in transmission, as shown in Fig. 1.

After passing through a phase modulating hologram, the electric field in the input plane has a modified wavefront

$$E^{in}(\stackrel{}{r})=E_0(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }^{in}(\stackrel{}{r})],$$
(8)
where $`\mathrm{\Phi }^{in}(\stackrel{}{r})`$ is the imposed phase profile. Calculating the phase hologram, $`\mathrm{\Phi }^{in}(\stackrel{}{r})`$, needed to project a desired pattern of traps is not particularly straightforward, as a simple example demonstrates.

In a typical application of holographic optical tweezer arrays, the undiffracted beam, $`E_0(\stackrel{}{r})`$, projects a single optical tweezer into the center of the focal plane with output wavefront $`E_0^f(\stackrel{}{\rho })`$, and the goal is to create displaced copies of this tweezer in the focal plane. One possible wavefront describing an array of $`N`$ optical tweezers at positions $`\stackrel{}{\rho }_i`$ in the focal plane is a superposition of single (non-overlapping) tweezers

$$E^f(\stackrel{}{\rho })=\underset{i=1}{\overset{N}{}}\alpha _iE_0^f(\stackrel{}{\rho }\stackrel{}{\rho }_i),$$
(9)
where the normalization $`_{i=1}^N|\alpha _i|^2=1`$ conserves energy. $`E^f(\stackrel{}{\rho })`$ may be written as a convolution

$`E^f(\stackrel{}{\rho })`$ $`={\displaystyle d^2\rho ^{}E_0^f(\stackrel{}{\rho }^{})T(\stackrel{}{\rho }\stackrel{}{\rho }^{})}`$ (10)
$`E_0^fT(\stackrel{}{\rho })`$ (11)
of $`E_0^f(\stackrel{}{\rho })`$ with a lattice function

$$T(\stackrel{}{\rho })=\underset{i=1}{\overset{N}{}}\alpha _i\delta ^{(2)}(\stackrel{}{\rho }\stackrel{}{\rho }_i).$$
(12)
Equations (6) and (8) relate $`E^f(\stackrel{}{\rho })`$ to the associated input wavefront:

$$\begin{array}{c}E_0^{in}(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }^{in}(\stackrel{}{r})]=^1\left\{E_0^fT(\stackrel{}{\rho })\right\}\hfill \\ \hfill =\frac{2\pi f}{k}^1\left\{E_0^f(\stackrel{}{\rho })\right\}^1\left\{T(\stackrel{}{\rho })\right\},\end{array}$$
(13)
by the Fourier convolution theorem. The phase modulation needed to achieve the array of optical tweezers then follows from Eq. (8):

$$\mathrm{exp}[i\mathrm{\Phi }^{in}(\stackrel{}{r})]=\frac{2\pi f}{k}^1\left\{T(\stackrel{}{\rho })\right\},$$
(14)
independent of the form of the single tweezer.

The phases of the complex weights, $`\alpha _i`$, must be selected so that $`\mathrm{\Phi }^{in}(\stackrel{}{r})`$ is a real-valued function. Unfortunately, the resulting system of equations has no analytic solution. Still greater difficulties are encountered in designing more general systems of optical traps, including tweezers which trap out of the focal plane or mixed arrays of conventional and vortex tweezers. Rather than deriving solutions for particular tweezer configurations, we have developed more general numerical methods which we apply in the following Sections to creating planar arrays optical tweezers.

## IV Adaptive-Additive Algorithm

Our approach is based on the adaptive-additive (AA) algorithm of Soifer *et al.* , an iterative numerical technique which explores the space of degenerate phase profiles, $`\mathrm{\Phi }^f(\stackrel{}{\rho })`$, to find a phase modulation of the incident laser beam encoding any desired intensity profile in the focal plane. To facilitate calculation and fabrication, both the input and output planes are discretized into $`M\times M`$ square arrays of pixels. Optimal spatial resolution requires pixels in the focal plane to be one half-wavelength on a side, $`\delta ^f=\lambda /2`$. The number, $`M`$, of pixels on a side then depends on the desired dimensions of the trapping array. Lengths in the input and focal planes are related by Eqs. (4) and (6), so that the corresponding pixel size in the input plane is $`\delta ^{in}=\lambda f/(M\delta ^f)=2f/M`$. If $`\delta ^{in}`$ is inconveniently small, then L1 and L2 can be chosen so that a more amenable pixel size at B corresponds to $`\delta ^{in}`$ at B.

The AA algorithm, depicted in Fig. 3, starts with an arbitrary initial guess for $`\mathrm{\Phi }_1^{in}(\stackrel{}{r})`$ and an initial input wavefront $`E_1^{in}(\stackrel{}{r})=E_0(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }_1^{in}(\stackrel{}{r})]`$. The Fourier transform of this wavefront is the starting estimate for the output electric field: $`E_1^f(\stackrel{}{\rho })=\{E_1^{in}(\stackrel{}{r})\}=A_1^f(\stackrel{}{\rho })\mathrm{exp}[i\mathrm{\Phi }_1^f(\stackrel{}{\rho })]`$. The corresponding intensity in the output plane, $`I_1^f(\stackrel{}{\rho })=|A_1^f(\stackrel{}{\rho })|^2`$ is unlikely to be a good rendition of the desired intensity pattern, $`I^f(\stackrel{}{\rho })=|A^f(\stackrel{}{\rho })|^2`$. The error,

$$ϵ_1\frac{1}{M^2}\underset{i=1}{\overset{M^2}{}}[I^f(\stackrel{}{\rho }_i)I_1^f(\stackrel{}{\rho }_i)]^2,$$
(15)
is reduced by mixing a proportion, $`a`$, of the desired amplitude into the field in the focal plane:

$$\overline{E}_1^f(\stackrel{}{\rho })=[aA^f(\stackrel{}{\rho })+(1a)A_1^f(\stackrel{}{\rho })]\mathrm{exp}[i\mathrm{\Phi }_1^f(\stackrel{}{\rho })].$$
(16)
Inverse transforming $`\overline{E}_1^f(\stackrel{}{\rho })`$ yields the corresponding field in the input plane, $`\overline{E}_1^{in}(\stackrel{}{r})=\overline{A}_1^{in}(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }_2^{in}(\stackrel{}{r})]`$. At this point, the amplitude in the input plane no longer matches the actual laser profile, so we replace $`\overline{A}_1^{in}(\stackrel{}{r})`$ with $`A_0(\stackrel{}{r})`$. The result is an improved estimate for the input field: $`E_2^{in}(\stackrel{}{r})=E_0(\stackrel{}{r})\mathrm{exp}[i\mathrm{\Phi }_2^{in}(\stackrel{}{r})]`$. This completes one iteration of the AA algorithm. Subsequent iterations lead to monotonically improving estimates, $`\mathrm{\Phi }_n^{in}(\stackrel{}{r})`$, for the desired phase modulation . The cycle is repeated until the error, $`ϵ_n`$, in the $`n`$-th iteration converges to within an acceptable tolerance: $`(ϵ_nϵ_{n1})/ϵ_n<\chi `$.

The phase and amplitude fields are computed as arrays of double-precision numbers, and their Fourier transforms calculated with fast Fourier transform (FFT) routines. Starting from random input phases, $`\mathrm{\Phi }_1^{in}(\stackrel{}{r}_i)`$, uniformly distributed in the range 0 to $`2\pi `$, the AA algorithm typically requires eight iterations to converge within $`\chi =10^6`$ of an acceptably accurate local minimum of $`ϵ_n`$ using an intermediate value for the mixing parameter, $`a=0.5`$.

## V Practical Considerations

The AA algorithm generates phase profiles, $`\mathrm{\Phi }^{in}(\stackrel{}{r})`$, that vary continuously between 0 and $`2\pi `$. Actually creating a phase element with continuously varying phase delay is difficult; usually only a small number of discrete levels are available. Discretizing the output of the adaptive-additive algorithm necessarily introduces errors. These can be minimized by integrating the discretization step into the AA algorithm itself, as shown in Fig. 3, although this can lead to problems with convergence.

### A Binarization

The most straightforward phase modulators offer just two levels of phase delay, and are known as binary holograms. Beyond quantization errors and their attendant loss of efficiency, binarization also imposes inversion symmetry on the output wavefront, $`E^f(\stackrel{}{\rho })=E^f(\stackrel{}{\rho })`$, and so limits what patterns can be generated. This might not seem a problem for inversion-symmetric patterns, but interference between two sides of the pattern can lead to unsatisfactory results, as shown in Figs. 4(a) and (b). If, however, we anticipate the reflection and calculate a phase mask encoding only half of the array, we achieve much better results, as shown in Figs. 4(c) and (d). In practice, we repeat this calculation about twenty times and choose the binary hologram with the best performance.

### B Tiling

Because fast Fourier transforms yield periodic functions, all holograms calculated with the AA algorithm can be tiled smoothly. That is, they can serve as the unit cell for new holograms without introducing phase discontinuities at the unit cell boundaries. The result of such tiling is to increase the spacing between tweezers by an amount proportional to the number of tilings along each dimension without reducing the resolution or trapping ability of the individual tweezers. Fig. 5 shows successive tilings of a hologram that generates a $`3\times 3`$ square array of tweezers, each labelled by the number of unit cells tiled along each side of the hologram. The same number describes the relative spacings of the resulting tweezers.

We use this property to design holograms encoding tweezer arrays with large inter-tweezer spacings. Increasing an array’s lattice constant requires smaller features in the input plane. In order to resolve these small features, the hologram’s pixel size must be reduced. Since the width of the hologram is fixed, the number of pixels increases with the inverse square of the pixel size. Calculating such holograms can become computationally expensive. Instead of directly calculating the hologram for a desired lattice constant, therefore, we calculate the smaller hologram encoding the same pattern with a proportionately smaller lattice constant. This hologram can be tiled to create a hologram for the desired tweezer spacing. Tiling can be done either numerically or physically, via a step and repeat mask fabrication process.

## VI Fabrication

Phase profiles can be recorded in the surface topography of an optical element , or in controlled variations in a dielectric’s index of refraction . Liquid crystal displays also have been used as phase-modulating elements , and dynamically reconfigurable tweezer arrays have been demonstrated in principle , although not yet in practice. Some photorefractive elements such as those being explored as optical memory devices also can be reconfigured, but must be programmed optically. Few, if any, are available as commercial optical elements. Photorefractive holograms created with photographic techniques  promise the greatest flexibility for creating static tweezer arrays at very low cost, but do not appear to have advanced beyond the research stage.

Surface patterning takes advantage of well-established photolithographic techniques and can be implemented easily and inexpensively. We have taken this approach in creating our own holographic optical tweezer arrays. Fig. 6 shows the principle. Light propagates more slowly in a dielectric material than in air. When a wavefront first enters the material, it is uniformly slowed to a speed $`c/n`$, where $`c`$ is the speed of light in vacuum and $`n`$ is the material’s index of refraction. Parts of the wavefront emerging first from the textured surface propagate at speed $`c`$, while sections remaining in the material fall behind, picking up a phase delay proportional to the extra thickness of material. Consequently, the relative phase at $`\stackrel{}{r}`$ is proportional to the surface’s relief, $`d(\stackrel{}{r})`$:

$$\mathrm{\Phi }^{in}(\stackrel{}{r})=2\pi (n1)\frac{d(\stackrel{}{r})}{\lambda }.$$
(17)
A similar principle applies when imposing a pattern of phase delays through the relief on a reflective surface, but with the factor $`n1`$ replaced by 2.

The pattern of hills and valleys needed to create a desired phase profile can be formed in photoelastic polymer gels. Such materials provide the recording medium for commercial holographic printers. These are not so common as photolithographic facilities for surface etching, however, so we digress in the next Section to describe the details of our fabrication process.

### A Reactive Ion Etching of Fused Silica

We employ reactive ion etching to create binary holograms in 1-mm thick substrates of polished fused silica, a transparent medium with an index of refraction of $`n=1.456`$ at the wavelength of our trapping laser, $`\lambda =532`$ nm. Obtaining a phase shift of $`\pi `$ radians requires a feature depth of $`\lambda /(2n2)=583`$ nm. The fabrication process has three main steps: creating a photomask, transferring this pattern to an etch mask covering the silica, and then etching to a precise depth, Fig. 7.

High contrast, high resolution film can be used to create masks for many of the trapping patterns we have investigated. We begin by laser printing the calculated phase profile as a binary image, with black pixels representing a relative phase shift of $`\pi `$ radians, and white representing 0 radians. This image is photoreduced to the actual dimensions of the hologram. Each of our holograms covers a square whose width, $`2f=3.24`$ mm, matches the laser beam’s diameter at B in Fig. 1. Holograms involving finer linewidths were vectorized before processing with commercial mask writers at the National Nanofabrication Facility.

We next create an etch mask on the surface of the fused-silica substrate. First the surface is protected with a 25 nm layer of chromium and a 1.76 $`\mu `$m layer of positive photoresist, Fig. 7(a). The photomask is placed in contact with the photoresist, and the entire sample is exposed to UV radiation, Fig. 7(b). The photomask is removed and exposed regions of the photoresist are dissolved away, revealing parts of the chromium layer, Fig. 7(c). Finally, the exposed chromium is removed with an acid wash, exposing the sections of silica to be etched, Fig. 7(d).

Unprotected regions of the silica are susceptible to attack by fluoride ions. Reactive ion etching provides a controlled exposure to ions generated by RF dissociation of a mixture of oxygen and carbon tetrafluoride. These reactive ions rapidly oxidize the organic photoresist, but are halted by the layer of metallic chromium. The unprotected regions of the silica surface continue to be removed at a rate of about 0.5 nm/sec, Fig. 7(e), until the etched regions reach the desired depth. As the final step, we remove the remaining chromium to reveal a precisely textured fused-silica surface, Fig. 7(f).

The etching process could be repeated with different photomasks to produce a more nuanced pattern . $`N`$ such steps would yield $`2^N`$ gradations of phase delay. Each step, however, would require planarizing and polishing the previously etched pattern, recoating the surface, and precisely aligning the new photomask over the existing pattern before etching. Not only is this is time consuming, it is not necessary for many applications.

### B Tolerances

Regardless of the fabrication method, any practical phase hologram will deviate from its design and these deviations will degrade its performance. We consider two principal fabrication defects: overall multiplicative error in the phase modulation due to mismatches between wavelength and etch depth, and random noise in the local phase shift due to roughness. To quantify these defects’ influence on hologram performance, we define the efficiency, $``$, to be the fraction of incident laser power projected into the planned tweezer pattern. For simplicity, we compare the intensity pattern in the focal plane when the actual hologram is illuminated by a uniform plane wave, $`\stackrel{~}{I}^f(\stackrel{}{\rho })`$, to the ideal intensity pattern in the focal plane, $`I_0^f(\stackrel{}{\rho })=(2\pi f/k)^2T^2(\stackrel{}{\rho })`$. The corresponding efficiency,

$$\frac{_{i=1}^{M^2}T^2(\stackrel{}{\rho }_i)\stackrel{~}{I}^f(\stackrel{}{\rho }_i)}{_{i=1}^{M^2}T^2(\stackrel{}{\rho }_i)},$$
(18)
is a less stringent measure of the agreement between the ideal and actual holograms than the error, $`ϵ_n`$, since it is possible to have $`=1`$ when $`ϵ_n>0`$, but $`ϵ_n=0`$ implies $`=1`$.

To give a feel for the results obtained with our methods, we calculate the efficiency of four standard holograms as a function of the magnitude of the fabrication defects. The four standard holograms are continuous and binary versions of holograms encoding $`4\times 4`$ and $`20\times 20`$ square tweezer arrays, each with the same lattice constant. We calculated all four holograms twenty times, and selected the most efficient hologram from each group to use in the the efficiency studies.

The phase modulation created by an etched hologram is proportional to the etch depth, Eq. (17). If the etch rate is not precisely controlled, or if the hologram is illuminated with light of the wrong wavelength, the actual phase profile, $`\stackrel{~}{\mathrm{\Phi }}^{in}(\stackrel{}{r})`$, will differ from the design $`\mathrm{\Phi }^{in}(\stackrel{}{r})`$ by a scale factor, $`\stackrel{~}{\mathrm{\Phi }}^{in}(\stackrel{}{r})=\alpha \mathrm{\Phi }^{in}(\stackrel{}{r})`$.

As $`\alpha `$ departs from unity, most of the laser light not contributing to the tweezer array is focused at the central undiffracted spot. Fig. 8 shows the efficiency of the four standard holograms as a function of $`\alpha `$. Even the continuous holograms with $`\alpha =1`$ are not perfectly efficient because the AA algorithm rarely identifies a globally ideal phase modulation. Binary holograms are still less efficient, with ideal efficiencies near 80%. Reassuringly, Fig. 8 suggests that a hologram’s efficiency does not depend strongly on precisely matching etch depth to the light’s wavelength.

Reactive ion-etching creates a rough surface, whose asperities add random fluctuations to the phase profile. We measured the surface topography of our fused silica wafers after etching and found a Gaussian distribution of etch depths, with a standard deviation of 60 nm or $`\pi /10`$ radians at 532 nm illumination. This roughness is laterally uncorrelated down to length scales of less than 280 nm.

We modeled such roughness by adding uncorrelated Gaussian noise to the calculated phase profiles,

$$\stackrel{~}{\mathrm{\Phi }}^{in}(\stackrel{}{r})=\mathrm{\Phi }^{in}(\stackrel{}{r})+\eta (\stackrel{}{r}),$$
(19)
where the noise’s probability distribution is given by

$$\rho (\eta )=\frac{1}{\sqrt{2\pi \sigma _\mathrm{\Phi }^2}}\mathrm{exp}\left(\frac{\eta ^2}{2\sigma _\mathrm{\Phi }^2}\right).$$
(20)
Fig. 9 shows how the efficiency of the four standard holograms decreases with increasing surface roughness.

Combining Eqs. (4) and (19) yields an expression for the electric field profile in the focal plane given a particular manifestation of the noise profile in the input plane,

$$\stackrel{~}{E}^f(\stackrel{}{\rho })=\frac{k}{2\pi f}d^2r\mathrm{exp}\left(i\frac{k}{f}\stackrel{}{r}\stackrel{}{\rho }+i\mathrm{\Phi }^{in}(\stackrel{}{r})+i\eta (\stackrel{}{r})\right).$$
(21)
Averaging over all possible phase profiles yields

$$\stackrel{~}{E}^f(\stackrel{}{\rho })=\mathrm{exp}\left(\frac{1}{2}\sigma _\mathrm{\Phi }^2\right)E^f(\stackrel{}{\rho }),$$
(22)
so that

$$(\sigma _\mathrm{\Phi })=(0)\mathrm{exp}\left(\sigma _\mathrm{\Phi }^2\right).$$
(23)
This result agrees well with numerically calculated efficiencies, as can be seen in Fig. 9. Substituting the measured $`\sigma _\mathrm{\Phi }^2`$ for our etched binary holograms, we estimate that roughness diminishes their efficiencies by a further 10% to roughly 70%.

## VII Further Considerations

Using the techniques described above, we have created triangular and square tweezer arrays which trap up to 400 particles at once. Still larger arrays and less regular arrangements are certainly feasible. Even static holograms permit some degree of reconfigurability. Rotating a hologram about its optical axis rotates the pattern of tweezers in the plane. Tilting it changes the aspect ratio. Individual traps can be turned off by blocking their beams in the plane conjugate to the object plane, labelled OP in Fig. 1. Such spatial filtering also can be useful for eliminating stray laser light, and to block out any undiffracted portion of the input beam. Replacing lenses L1 and L2 with zoom lenses should permit a degree of continuous scaling of the lattice constant.

The methods described in the previous Sections are appropriate for projecting arrays of identical tweezers in the plane, where each tweezer shares the properties of a single tweezer formed by the unmodulated input beam. Shaping the wavefronts of the individual beams, for example to embed some optical vortices in an array of conventional optical tweezers, requires a straightforward elaboration of the AA algorithm . Creating three dimensional arrays, on the other hand, requires more sophisticated calculations to avoid undesirable interference effects, and will be discussed elsewhere.

## VIII Acknowledgements

This work was funded by the National Science Foundation through grant DMR-9730189, by the MRSEC program of the NSF through grant DMR-9888595 and by the David and Lucile Packard Foundation.


format: text

text:
# An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions

## Abstract

The energy and baryon densities in heavy ion collisions are estimated by analytically solving a 1+1 dimensional hydrodynamical model with source terms. Particularly, a competition between the energy and baryon sources and the expansion of the system is discussed in detail.

PACS: 24.10.Nz, 25.75.-q, 12.38.Mh

Theoretical evidence from lattice simulations of Quantum chromodynamics (QCD) for a phase transition from a hadron gas to a quark-gluon plasma at high temperatures has prompted experimental efforts to create this new phase of matter in the laboratory during the early stages of relativistic heavy-ion collisions. Besides the deconfinement phase transition, the medium effects in hot and dense nuclear matter produced in the collisions change the characteristics of strong interaction at hadronic level. For instance, the effects lead to the transition from the chiral symmetry breaking phase to the phase in which it is restored. Since 1986, oxygen, sulfur, silicon and other beams have been successfully accelerated to 14.5 A GeV at the BNL AGS, and to 200 A GeV at the CERN SPS. In the past few years heavy-ion collisions are also investigated at the GSI SIS with relative low energies of 0.8-1.8 A GeV to study the in-medium effects related to high densities. From the above experiments, people have really seen some new phenomena such as the $`J/\psi `$-suppression, the dilepton distribution shift and the $`\frac{K^{}}{K^+}`$ enhancement. It is, however, still difficult to say that definite signals of formation of quark-gluon plasma in the collisions are found. This vague situation is expected to be changed in the near future, since the BNL RHIC will start to run soon, providing center of mass energies of up to 100 + 100 A GeV and beams as heavy as uranium. To further study the high density effects at hadronic level, a new generation of Cooling Storage Ring(CSR) with incident energies of 0.4-0.9 A GeV is going to be constructed at LanZhou, China.

The key point for producing quark-gluon plasma or discussing in-medium effects in heavy ion collisions is whether the energy density or baryon density achieved in the collisions is high enough. From a general consideration, the energy density required for producing the quark-gluon plasma must be several times higher than the energy density of normal nuclear matter, and the baryon density for a visible change of chiral properties in nuclear matter should be at least the order of the normal nuclear density. Can such a high energy density or baryon density be achieved at SPS and RHIC or at AGS, SIS and CSR? The answer to this question depends strongly on the degree of nuclear stopping power. The experimental data from AGS demonstrated quantitatively that around and below the AGS energy any incident nucleus can be fully stopped in a heavy target. In contrast with this, the data from SPS show that the incident energy is not fully deposited in the target, the colliding nuclei are partially transparent. Since the degree of stopping power is variable with the incident energy, it is believed that a baryon-rich matter may be created either in the central region at lower energies or in the fragmentation region at higher energies.

Assuming full stopping, the maximum energy density and baryon density in the center of mass frame for a symmetrical $`AA`$ collision can be simply evaluated to be

$`ϵ_{max}=2n_Nm_N\mathrm{cosh}^2(y_c),`$
$`n_{max}=2n_N\mathrm{cosh}(y_c),`$ (1)
where $`n_N`$ is the density of nucleons in normal nuclei, $`m_N`$ the mass of the nucleon, and $`y_c`$ the rapidity of the colliding nuclei in the center of mass frame.

For the collisions with partial transparency, Bjorken’s formula gives

$$ϵ_{max}=\frac{1}{\pi R_A^2\tau _0}\left(\frac{dE_T}{dy}\right)_{max},$$
(2)
where $`\tau _0`$ is the formation time of the secondary particles, $`R_A`$ the radius of the colliding nuclei, and $`\left(\frac{dE_T}{dy}\right)_{max}`$ the transverse energy distribution in the central region.

Hydrodynamical models are widely used to describe the macroscopic and collective effects in heavy ion collisions, assuming thermal and chemical equilibrium through a large number of nucleon-nucleon and secondary collisions. Landau was the first to treat high energy collisions with fluid dynamics, and Bjorken considered partial transparency in relativistic collisions. Assuming the colliding system is boost invariant, the maximum energy density reached at time $`\tau _0`$ decreases monotonously as

$$ϵ(\tau )=ϵ_{max}\left(\frac{\tau _0}{\tau }\right)^{1+v_s^2}$$
(3)
in the central rapidity region, where $`v_s`$ is the sound velocity in the fluid.

In Landau and Bjorken’s discussions there is no consideration of the initial collision process. Two nuclei are assumed to collide at a unique space-time point $`(x,t)=(0,0)`$. This means that all nucleons in the projectile collide with the target at the same time $`\tau =0`$, and the secondary particles from different collisions are formed at the same time $`\tau =\tau _0`$. Therefore, the maximum energy density $`ϵ_{max}`$ in (An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions) or (2) is just a simple superposition of the $`NA`$ collisions. After time $`\tau _0`$, the energy density drops down in terms of (3) due to the expansion of the system.

It is clear that any nucleus-nucleus ($`AA`$) collision with a finite rapidity $`y_c`$ in the center of mass frame takes a finite collision time

$$d_t=\frac{2(R_Ar_N)}{\mathrm{sinh}y_c}0,$$
(4)
long or short according to the incident energy and the size of the colliding nuclei. Here $`r_N`$ is the nucleon radius. Considering a central $`{}_{}{}^{197}Au^{197}Au`$ collision at incident energies $`E_{lab}=0.414.5`$ A GeV, namely from CSR to AGS energies, the collision time $`4fmd_t23fm`$ is comparable with or even larger than the total time of the fluid description of the system. However, when we consider the same colliding nuclei but at RHIC energy, the collision time is only $`0.1fm`$, the neglect of its influence on the fluid evolution is therefore a good approximation in high energy limit. In this paper, we will concentrate our attention on the dependence of the energy and baryon densities on the collision time for heavy ion collisions with full nuclear stopping.

The boost invariance in the longitudinal direction used to derive the scaling solution (3) is a good assumption at early stages of high-energy nuclear collisions when the transverse expansion can be neglected. As the transverse flow develops with increasing time, the fluid cools down more and more rapidly. At the late time in the evolution of the reaction, the fluid will expand radially outward, and one can use boost invariance in the radial direction instead of in the longitudinal direction. In fact, the radial scaling solution is widely used in analytic study of relativistic heavy-ion collisions, for instance, the investigation of disoriented chiral condensate. Compared with the high-energy collisions, the kinematics dominated longitudinal expansion in low-energy collisions is much weaker, but the transverse expansion does not change a lot. For example, when the colliding energy increases from about $`20`$ A GeV at SPS to $`200`$ A GeV at RHIC, the transverse flow velocity changes from $`0.55c`$ to $`0.65c`$. On the other hand, what we are interested in here is the initial collision period. At high energies, this period is very short, the longitudinal scaling solution is a good approximation in the period. At low energies, the collision time $`d_t`$ has the same order as the total evolution time of the fluid. For most time of the collision period, the transverse flow is comparable with the longitudinal flow. Therefore, for analytically estimating the energy density and baryon number density in low-energy collisions the radial scaling solution seems better than the longitudinal scaling solution. In this case the time evolution of the energy density (3) is replaced by

$$ϵ(\tau )=ϵ_{max}\left(\frac{\tau _0}{\tau }\right)^{3(1+v_s^2)},$$
(5)
the factor $`3`$ in the exponent arises from the 3-dimensional homogeneous expansion.

In principle, the expansion of the system starts at the initial time $`\tau _0`$. However, in the very early stage the expansion is weaker than the homogeneous evolution (5). For the sake of simplicity we still use the homogeneous mechanism as an effective description of the whole expansion process, but the beginning time of the expansion is shifted from $`\tau _0`$ to $`\tau _{exp}`$. Since the collision time $`d_t`$ is large enough, $`\tau _{exp}`$ can be restricted in the region $`\tau _0\tau _{exp}\tau _0+d_t`$.

Let us first analyze qualitatively the behavior of the densities when a non-zero collision time is taken into account. We still use $`NA`$ collisions as input, but the collision time of each $`NA`$ collision is different. This implies a different time in the formation of secondary particles. Those produced by the first $`NA`$ collision are formed at time $`\tau _0`$, those by the last one at $`\tau _0+d_t`$, and others between $`\tau _0`$ and $`\tau _0+d_t`$. As there is a continuous input of secondaries and no expansion before $`\tau _{exp}`$, the energy and baryon densities of the system will increase linearly with time in the range $`\tau _0\tau \tau _{exp}`$. When the effective expansion begins, there is a competition between the input from the left $`NA`$ collisions and the expansion of the system. While the former makes the densities increase, the latter attenuates the fluid strongly. Therefore, the evolution of the densities in the time interval $`\tau _{exp}\tau \tau _0+d_t`$ depends on the input and expansion rates. After the time $`\tau _0+d_t`$ when the input from the $`NA`$ collisions ends, the densities drop down monotonously like Eq.(5).

In the absence of viscosity and thermal conductivity, the behavior of a relativistic fluid is controlled by the energy-momentum and baryon number conservation laws,

$`_\mu T^{\mu \nu }=\mathrm{\Sigma }^\nu ,`$
$`_\mu n^\mu =\sigma ,\mu ,\nu =0,1,2,3`$ (6)
where $`T^{\mu \nu }=Pg^{\mu \nu }+(ϵ+P)u^\mu u^\nu `$ is the energy-momentum tensor, $`n^\mu =nu^\mu `$ the baryon flow, $`P=P(ϵ)`$ the pressure, $`ϵ`$ the local energy density and $`u^\mu `$ the four-velocity of the fluid. The source terms $`\mathrm{\Sigma }^\mu `$ and $`\sigma `$ indicate the continuous deposition of energy, momentum and baryon number in the fluid due to the $`NA`$ collisions. In ultra-relativistic collisions, there is almost no net baryon number in the central region, one discusses the energy density only. For collisions with full nuclear stopping, we must consider simultaneously the energy and baryon densities characterized by (An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions).

Using two Lorentz covariant coordinate, the proper time $`\tau =\sqrt{t^2r^2}`$ and the space-time rapidity $`y=\frac{1}{2}\mathrm{ln}\frac{t+r}{tr}`$ instead of $`r`$ and $`t`$, the above hydrodynamical equations are greatly simplified as

$`_\tau ϵ+3{\displaystyle \frac{P+ϵ}{\tau }}=\mathrm{\Sigma },`$
$`_\tau n+3{\displaystyle \frac{n}{\tau }}=\sigma ,`$ (7)
where the source terms on the right side mean the energy density and baryon density deposited in the fluid per unit time. Since the total energy input and the total baryon number input during the collision time $`d_t`$ are $`ϵ_{max}(A^{1/3}1)/A^{1/3}`$ and $`n_{max}(A^{1/3}1)/A^{1/3}`$, where $`ϵ_{max}`$ and $`n_{max}`$ are shown in (An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions), the source terms are given by

$`\mathrm{\Sigma }={\displaystyle \frac{ϵ_{max}(A^{1/3}1)}{A^{1/3}d_t}},`$
$`\sigma ={\displaystyle \frac{n_{max}(A^{1/3}1)}{A^{1/3}d_t}}.`$ (8)
Taking into account the initial condition

$`ϵ(\tau _0)={\displaystyle \frac{ϵ_{max}}{A^{1/3}}},`$
$`n(\tau _0)={\displaystyle \frac{n_{max}}{A^{1/3}}},`$ (9)
the solution of the differential equations (An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions) is

$`ϵ(\tau )=\{\begin{array}{cc}ϵ(\tau _0)+(\tau \tau _0)\mathrm{\Sigma }\hfill & \tau _0\tau \tau _{exp},\hfill \\ ϵ(\tau _{exp})\left(\frac{\tau _{exp}}{\tau }\right)^{3(1+v_s^2)}+\frac{\mathrm{\Sigma }}{3(1+v_s^2)+1}\tau \left(1\left(\frac{\tau _{exp}}{\tau }\right)^{3(1+v_s^2)+1}\right)\hfill & \tau _{exp}\tau \tau _0+d_t,\hfill \\ ϵ(\tau _0+d_t)\left(\frac{\tau _0+d_t}{\tau }\right)^{3(1+v_s^2)}\hfill & \tau _0+d_t\tau ,\hfill \end{array}`$ (13)
$`n(\tau )=\{\begin{array}{cc}n(\tau _0)+(\tau \tau _0)\sigma \hfill & \tau _0\tau \tau _{exp},\hfill \\ n(\tau _{exp})\left(\frac{\tau _{exp}}{\tau }\right)^3+\frac{\sigma }{4}\tau \left(1\left(\frac{\tau _{exp}}{\tau }\right)^4\right)\hfill & \tau _{exp}\tau \tau _0+d_t,\hfill \\ n(\tau _0+d_t)\left(\frac{\tau _0+d_t}{\tau }\right)^3\hfill & \tau _0+d_t\tau ,\hfill \end{array}`$ (17)
where the sound velocity $`v_s`$ comes from the equation of state $`P=v_s^2ϵ`$. Before the effective expansion begins, the energy density increases linearly. In the time period $`\tau _{exp}\tau \tau _0+d_t`$, the energy density contains two terms. The first one is the time evolution of the density at $`\tau _{exp}`$, and the second term reflects the contribution of the source term. After $`\tau _0+d_t`$ there is no more energy input, the second term disappears and the density decreases with time. The behavior of the baryon density is similar to that of the energy density. Obviously the heavier the colliding nuclei are or the lower the incident energy is, the longer is the collision time $`d_t`$ and the more important is the energy and baryon diffusion.

Let us now turn to numerical results. As we have emphasized above, the effect of finite collision time on the density evolution is especially important for collisions with full stopping. In our numerical calculations below we take $`{}_{}{}^{197}Au^{197}Au`$ collisions at AGS, SIS and CSR energies, namely from $`0.4`$ to $`14.5`$ A GeV. At these low energies, while it seems not possible to create quark-gluon plasma in the collisions, the in-medium effects at hadronic level are still visible if the energy density or baryon density is larger than the relevant density of normal nuclear matter. For instance, from a model independent calculation the broken chiral symmetry in the vacuum $`(T=n=0)`$ is partially restored even in the normal nuclear matter. The energy density and baryon density that are shown in the following are scaled, respectively, by the normal energy density $`ϵ_N=m_Nn_N`$ and normal baryon density $`n_N=1/\left(\frac{4}{3}\pi r_N^3\right)`$, where the nucleon mass $`m_N`$ and the nucleon radius $`r_N`$ are taken to be $`0.94`$ GeV and $`1.1`$ fm. For a perfect fluid, the equation of state can be written as $`P=\frac{1}{3}ϵ`$. In real fluid with finite hadron mass, viscosity and thermal conductivity, one has $`P=v_s^2ϵ`$ with $`v_s^2<\frac{1}{3}`$. In our calculation we take $`v_s^2=\frac{1}{5}`$.

In Fig.1a, the energy density, scaled by $`ϵ_N`$, is shown as a function of the proper time $`\tau `$ for four values of the expansion time $`\tau _{exp}`$ at the incident energy $`E_{lab}=1.8`$ A GeV. For $`\tau _{exp}=\tau _0`$, the effective expansion and the energy input start simultaneously. The energy density drops down in the very beginning when the energy input rate is less than the expansion rate, and then goes up smoothly when the energy deposited in the fluid is high enough. The maximum energy density $`ϵ(\tau _0+d_t)/ϵ_N=0.76`$ is larger than the initial value $`ϵ(\tau _0)/ϵ_N=0.67`$. When the fluid evolves beyond the source region, the energy density decreases monotonously. As the effective expansion starts later than the formation of the fluid, $`\tau _{exp}>\tau _0`$, the energy density increases linearly before the expansion. The increase or decrease of the energy density in the region $`\tau _{exp}\tau \tau _0+d_t`$ depends on the competition between the amount of the energy deposited in the fluid and the expansion rate. When $`\tau _{exp}`$ is larger than $`1.3\tau _0`$, the position of the maximum energy density moves from $`\tau _0+d_t`$ to $`\tau _{exp}`$. For $`\tau _{exp}\tau _0+d_t`$, there is no expansion in the source region, the maximum energy density $`ϵ(\tau _0+d_t)`$ achieves its limit value (An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions).

The time evolution of the baryon density, scaled by $`n_N`$, is shown in Fig.1b. The structure is very similar to that of the energy density.

In Fig.2, we show the maximum energy density and baryon density, scaled by $`ϵ_N`$ and $`n_N`$, as functions of incident energy $`E_{lab}`$ for the four values of the expansion time $`\tau _{exp}`$. In the case of no expansion in the source region , the maximum energy density and baryon density given by (An Analytic Solution of Hydrodynamic Equations with Source Terms in Heavy Ion Collisions) increase very rapidly and reach $`17.5`$ and $`6`$ at $`E_{lab}=14.5`$ A GeV. When the expansion in the source region is included, the maximum densities are strongly suppressed by the expansion.

Obviously the energy and baryon densities in heavy ion collisions can not reach the values estimated in the case of no expansion in the source region. The attenuation of the fluid in the initial collision period must be considered. As is well known, the in-medium effects in low energy collisions are mainly related to the baryon number density. When the in-medium effects become visible, the baryon density should be at least larger than the normal nuclear density, $`n_{max}/n_N1`$. The experimental data of heavy ion collisions at SIS with incident energies $`E_{lab}=0.81.8`$ A GeV, such as the $`K^{}/K^+`$ enhancement, do show the importance of the in-medium effects. This means that the baryon density at SIS energies is already beyond the normal nuclear density. From the comparison of this analysis with our calculations shown in Figs.1 and 2, we can choose approximately $`\tau _0+d_t/2`$ as the beginning time of the effective expansion for the estimation of the energy and baryon densities in low energy collisions.

In summary, we have investigated the influence of non-zero collision time on the evolution of energy and baryon densities in heavy ion collisions. When the size of the colliding nuclei is large enough and the incident energy is not very high, for instance the collisions at AGS, SIS and CSR, the initial collision time is long, and the colliding energy is deposited in the fluid gradually instead of an instantaneous input in ultra-relativistic collisions at RHIC and LHC. In the collision period, there exists a competition between the energy and baryon deposition and the expansion. The densities inside and outside the source region behavior very differently. They are strongly suppressed by the expansion in the source region, and the time to reach the maxima is also remarkably delayed.

The work was supported in part by the NSFC under grant numbers 19845001 and 19925519, and by the Major State Basic Research Development Program under contract number G2000077407.

FIGURE CAPTIONS

Fig.1
The time evolution of the energy density (Fig.1a) and baryon density (Fig.1b), scaled by $`ϵ_N`$ and $`n_N`$, in $`{}_{}{}^{197}Au^{197}Au`$ collisions at incident energy $`E_{lab}=1.8`$ A GeV for four values of the expansion time $`\tau _{exp}`$.

Fig.2
The maximum energy density (Fig.2a) and baryon density (Fig.2b), scaled by $`ϵ_N`$ and $`n_N`$, as functions of incident energy $`E_{lab}`$ for four values of the expansion time $`\tau _{exp}`$.