text:
# LiBC by Polarized Raman Spectroscopy: Evidence for Lower Crystal Symmetry?

## Abstract

The paper presents polarized Raman scattering study on a few-micron-size crystallite of LiBC with natural faces. The experiment on as grown sample has revealed a four lattice modes with frequencies at 1276 cm-1, 830 cm-1, 546 cm-1 and 170 cm-1, respectively. The number of observed Raman lines and their selection rules are incompatible with the assumed $`\mathrm{D}_{6\mathrm{h}}`$ symmetry. The modes at 1276 cm-1 and 170 cm-1 correspond to the expected Raman active modes. In contrast with the superconducting compound MgB2, the B-C bond stretching mode (at 1276 cm-1) has rather small damping. The two ”forbidden” modes (at 830 cm-1 and 546 cm-1) disappeared after subsequent thermal treatment.

Recently there were many attempts to further increase the remarkably high phase transition temperature toward the superconducting state in MgB2nature1  by various chemical substitutions and doping, but with very little progress.Ros02  A new hope was brought by theoretical proposition that in a related isoelectronic compound, LiBC, significantly higher phase transition temperatures can be attained by hole doping produced by Li non-stoichiometry or by field injection techniqueRos02 . Since the spectroscopic information on the phonon spectra in these electron-phonon coupling type superconductorsBCS ; Kortus ; Bela ; An ; BOH ; Yil  is obviously important, we have undertaken a basic polarized Raman scattering characterization of both MgB2my  and LiBC. The previous polarized Raman study on MgB2my  have proven that the unusually broad spectral feature around 600 cm-1 corresponds indeed to the $`\mathrm{E}_{2\mathrm{g}}`$ zone center mode. The present polarized Raman study complements the recent results on LiBC phonons by infrared spectroscopyArtem ; Bha  and ab-initio calculationsAn02 .

According to Ref.Wor95 , the LiBC crystals belongs to the $`\mathrm{D}_{6\mathrm{h}}^4`$ ($`\mathrm{P6}_3/\mathrm{mmc}`$) space group symmetry, with Li, B and C atoms in 2a, 2c and 2d Wyckoff positions, respectively. Such structure is thus similar to MgB2, except for the replacement of Mg by Li and by replacement of B by C at every second position along in-plane covalent bonds as well as along the hexagonal axis, what leads to doubled unit cell along the hexagonal axis. The factor group analysis at the $`\mathrm{\Gamma }`$-point yields $`2\mathrm{A}_{2\mathrm{u}}+2\mathrm{B}_{1\mathrm{g}}+\mathrm{B}_{2\mathrm{u}}`$ and $`2\mathrm{E}_{1\mathrm{u}}+2\mathrm{E}_{2\mathrm{g}}+\mathrm{E}_{2\mathrm{u}}`$ zone center optic modes with atomic displacements parallel and perpendicular to the hexagonal axis, respectively. Because of the unit cell doubling, the internal modes of B-C planes form Davydov-like doublets: one pair of B-C bond stretching modes ($`1\mathrm{E}_{1\mathrm{u}}+1\mathrm{E}_{2\mathrm{g}}`$) and a pair of ring puckering modes ($`1\mathrm{A}_{2\mathrm{u}}+1\mathrm{B}_{1\mathrm{g}}`$). The remaining six modes can be understood as ”external” modes including a pair of purely Li modes ($`1\mathrm{E}_{2\mathrm{u}}+1\mathrm{B}_{1\mathrm{g}}`$) and four other modes in which the B-C planes vibrate without deformation. Among all these modes, only the $`\mathrm{E}_{2\mathrm{g}}`$ representation is Raman active, so that only two first order Raman modes are to be expected in LiBC.

The experiments were carried out at room temperature, using a Renishaw Raman microscope with 514.5 nm (2.41 eV) argon laser excitation. The instrument allows measurements of polarized Raman spectra in back scattering from a spot size down to 1-2 microns in diameter. To minimize heating of the sample in the laser focus, the laser power was kept below 1mW. We have used the same, about 10 micron thick platelet sample, as in the recent IR reflectivity studyArtem . A regular hexagonal flat terrace with a diameter of about 10 microns (see Fig. 1 of the Ref. Artem, ) in the corner of the sample confirmsJung  that the hexagonal axis is perpendicular to the platelet sample.

The backscattering cross-polarized spectrum measured directly from this small hexagonal face is shown in the lower curve of the Fig. 1. As expected, there are two sharp Raman lines. Obviously, the higher frequency mode (near 1176 cm-1) corresponds to the $`\mathrm{E}_{2\mathrm{g}}`$ bond stretching mode (neighbors in B-C planes vibrate in anti-phase, neighbors along hexagonal axis vibrate in-phase, Li does not participate). It has practically the same frequency as the $`\mathrm{E}_{1\mathrm{u}}`$ bond stretching mode ( neighbors in B-C planes and neighbors along hexagonal axis vibrate in anti-phase), seen clearly in the IR reflectivity spectraArtem . This shows that the Davydov-like splitting is rather small, in agreement with the ab-initio calculationsAn02  predicting the two zone center bond stretching modes atBha  1185 cm-1 and 1194 cm-1. Let us note that unlike in MgB2, both these bond stretching modes in LiBC have rather small damping. Therefore, the electron-phonon interaction properties of LiBC are probably closer to AlB2.Pos02 ; Ren02

The eigenvector of the second expected $`\mathrm{E}_{2\mathrm{g}}`$ mode represents a shear mode of rigid B-C planes (all atoms in B-C plane move in the same sense, the next layer vibrates in anti-phase, Li does not participate). Such mode obviously derives from the transverse acoustic mode at the A-point of the parent MgB2 structure, so that its frequency can be estimated from the velocity of the transverse acoustic modes propagating along the hexagonal axis. The frequency of the corresponding A-point mode in MgB2 and AlB2 can be read from Fig. 3 of Ref. Ren02,  (200 cm-1 and 160 cm-1 respectively). Indeed, this estimation nicely agrees with the frequency of the other mode (near 170 cm-1) observed in our cross-polarized spectra shown in Fig.1.

The upper curve in the Fig. 1 was taken from the exactly same spot but with parallel polarization. In agreement with the expected symmetry properties of Raman tensors for the $`\mathrm{E}_{2\mathrm{g}}`$ modes, the intensity of the two above discussed modes does not change significantly. At the same time, several new lines have appeared in the spectrum. Somewhat higher luminescence background and broad features d1, d2 and d3 corresponding to the pronounced bands of LiBC phonon density states (see Fig. 3 of Ref. An02, , Fig. 4 of Ref. Bha, ) are rather common for parallel-polarized scattering configurationmy . However, the presence of the two additional sharp lines near 546 cm-1 and 830 cm-1 is in a flagrant disagreement with the above symmetry analysis, as if the actual factor group symmetry were lower than that of the expected $`\mathrm{D}_{6\mathrm{h}}`$ group.

There is no doubt that the 830 cm-1 mode is a ”ring-puckering” mode. It is too high frequency for an ”external” mode and the frequency region of bond stretching modes is limited by the d2 and d3 edges of the phonon density of states. Moreover, its frequency roughly coincides with ring-puckering phonon band calculated ab-initioAn02 .

There are certainly various mechanisms that may be invoked to explain the Raman activity of the ring-puckering mode in as-grown LiBC, but a decisive proof for the actual mechanism is in anyway beyond the scope of this study. Nevertheless, we would like to point out that a rather straightforward explanation can be obtained when assuming that the puckering of B-C planes is in fact present as a static distortion. In this case the puckering mode, which is ”frozen in” the parent $`\mathrm{P6}_3/\mathrm{mmc}`$ structure, becomes a totally symmetric mode of the distorted structure, and such a mode is always Raman active in parallel-polarized geometry. Since just four clear Raman modes have been observed in our spectra, we have searched possible candidates for frozen-in modes among the zone center modes only. As mentioned previously, there are only two such puckering modes ($`1\mathrm{A}_{2\mathrm{u}}+1\mathrm{B}_{1\mathrm{g}}`$). We can exclude the $`\mathrm{A}_{2\mathrm{u}}`$ mode because it would reduce the factor group symmetry towards $`\mathrm{C}_{6\mathrm{v}}`$ (within $`\mathrm{C}_{6\mathrm{v}}`$, the $`\mathrm{E}_{1\mathrm{u}}`$ should appear together with the two $`\mathrm{E}_{2\mathrm{g}}`$ modes as 3 $`\mathrm{E}_2`$ modes in the cross-polarized spectra, while only two modes were observed there). Thus, we are left with a single candidate, the $`\mathrm{B}_{1\mathrm{g}}`$ mode. The $`\mathrm{B}_{1\mathrm{g}}`$ zone center mode distortion leads to the $`\mathrm{D}_{3\mathrm{d}}^3`$ ($`\mathrm{P}\overline{3}\mathrm{m1}`$) structure, which is compatible with our experimental data, as we shall show in the following.

First of all, the factor group analysis for the hypothetic distorted structure yields $`3\mathrm{A}_{2\mathrm{u}}+2\mathrm{A}_{1\mathrm{g}}`$ and $`3\mathrm{E}_\mathrm{u}+2\mathrm{E}_\mathrm{g}`$ zone center optic modes with atomic displacements parallel and perpendicular to the trigonal axis, respectively. They are all either Raman ($`2\mathrm{A}_{1\mathrm{g}}`$ and $`2\mathrm{E}_\mathrm{g}`$) or infrared ($`3\mathrm{A}_{2\mathrm{u}}`$ and $`3\mathrm{E}_\mathrm{u}`$) active modes. To facilitate further discussion, from now on we distinguish the symmetry adapted eigenvectors by the symbols of correspondent irreducible representation in both $`\mathrm{D}_{6\mathrm{h}}`$ and $`\mathrm{D}_{3\mathrm{d}}`$ factor groups, the latter being in the parentheses.

The observation shown in the Fig. 1 is in agreement with the Raman selection rules for the $`\mathrm{D}_{3\mathrm{d}}`$ factor group. The pair of $`\mathrm{B}_{1\mathrm{g}}(\mathrm{A}_\mathrm{g})`$ modes should be observable in parallel-polarized but not in cross-polarized geometry, which indeed holds for the lines near 546 cm-1 and 830 cm-1. The requirement that the pair of $`\mathrm{E}_{2\mathrm{g}}(\mathrm{E}_\mathrm{g})`$ modes should be observable in both cross-polarized and parallel-polarized geometry with the same intensity holds both for $`\mathrm{D}_{3\mathrm{d}}`$ and $`\mathrm{D}_{6\mathrm{h}}`$) symmetry ( and is thus also fulfilled).

Measurements analogous to those shown in the Fig. 1 were reproduced at several places in the sample, both inside and outside the small hexagonal terrace, and for different orientations of incident light polarization axis, all with similar result. To obtain one more supplementary argument in favor of our conclusion, we have tried to measure the remaining polarization components from backscattering on the edges of the sample. Unfortunately, we have not found well developed flat surfaces perpendicular to the platelet, so that this measurement was rather difficult and the phonon propagation direction may not be strictly perpendicular to the c-axis. Nevertheless, the obtained results are again in agreement with the proposed assignment within the $`\mathrm{D}_{3\mathrm{d}}`$ factor group. The pair of $`\mathrm{E}_{2\mathrm{g}}(\mathrm{E}_\mathrm{g})`$ modes is present even when one of the polarization axes is parallel to the hexagonal (trigonal) axis (see ”zz” and ”zx” spectra shown in Fig. 2), which is expected within $`\mathrm{D}_{3\mathrm{d}}`$ but not allowed in $`\mathrm{D}_{6\mathrm{h}}`$. Similarly, the $`\mathrm{B}_{1\mathrm{g}}(\mathrm{A}_\mathrm{g})`$ modes, are missing in cross-polarized geometry, as required by $`\mathrm{D}_{3\mathrm{d}}`$ group Raman selection rules. Finally, no frequency shifts which could appear in the case of polar modes due to the LO-TO splitting were observed, again in agreement with our above conjecture that the frozen displacement is not the $`\mathrm{A}_{2\mathrm{u}}`$ polar mode.

Samples of the same origin were previously studied by infrared spectroscopy.Artem  Among the six remaining zone center optic modes, three have dipole moment along the trigonal axis ($`2\mathrm{A}_{2\mathrm{u}}(\mathrm{A}_{2\mathrm{u}})+1\mathrm{B}_{2\mathrm{u}}(\mathrm{A}_{2\mathrm{u}})`$) and three have dipole moments perpendicular to the trigonal axis ($`2\mathrm{E}_{1\mathrm{u}}(\mathrm{E}_\mathrm{u})+1\mathrm{E}_{2\mathrm{u}}(\mathrm{E}_\mathrm{u})`$). The latter triplet should contribute in the infrared reflectivity, measured from the large surface of the investigated sample. This is in agreement with the measurements of Ref. Artem, , which shows a clear B-C stretching mode near 1180 cm-1 (this has to be $`\mathrm{E}_{1\mathrm{u}}(\mathrm{E}_\mathrm{u})`$ mode) and a broad band with two well defined peaks (which could be tentatively assigned to the remaining $`\mathrm{E}_{1\mathrm{u}}(\mathrm{E}_\mathrm{u})+\mathrm{E}_{2\mathrm{u}}(\mathrm{E}_\mathrm{u})`$ modes).

After completing of the above described investigations, the sample was placed in a variable-temperature cell (LINKAM) in a hope to drive the system across a hypothetic phase transformation towards the prototype P63/mmc structure by thermal treatment. For this reason, the sample was first cooled down to about 80 K and then heated up in steps of about 50 K in order to see the evolution of the Raman spectrum. There was no clear indication of such a phase transformation up to 400 K. Unfortunately, in the vicinity of 450 K, a sudden drop-out of electricity stopped the experiment. Continuation of the experiment was possible two days later. This time the measurements were done while heating the sample from the room temperature up to 650 K. To our surprise, the two forbidden lines were never observed again.

In conclusion, we have observed scattering by the two $`\mathrm{E}_{2\mathrm{g}}`$ LiBC Raman active modes ( at 1276 cm-1, and 172 cm-1, respectively). The relatively high frequency and small damping of the stretching mode correlatesAn02  with the absence of the superconductivity in the stoichiometric LiBC. This indicates that the Raman scattering provides a very useful tool for easy testing of MgB2 related superconducting materials.

We have found that results of the polarized Raman measurements on the as-grown LiBC single crystal are incompatible with the assumed LiBC structure. Most unexpectedly, we have found two other modes at 830 cm-1 and 546 cm-1. The former forbidden mode is probably the $`\mathrm{B}_{1\mathrm{g}}`$ B-C puckering mode. The observed results can be explained assuming lower crystal symmetry, e.g $`\mathrm{P}\overline{3}\mathrm{m1}`$ space group symmetry induced by puckering displacement of B-C planes. This need not be the correct picture (ab-initio test calculations have not found such puckered structure stableprivate ), but it nevertheless gives some insight into the phenomenon.

The two forbidden Raman modes vanished after a subsequent thermal treatment. The influence of thermal treatment of LiBC based materials seems to be of essential importance for their superconductor properties, and the appearance (and disappearance) of forbidden Raman lines in LiBC samples may be one of the clue in understanding of the relevant microscopical processes. For example, B-C order stacking faults and/or Li non-stoichiometry allows a number of possible scenarios. Obviously, detailed understanding calls for further investigations. It seems that the sample is sensitive to heating to temperatures exceeding 400 K. We plan to perform Raman investigations of the thermal treatment in LiBC on new samples in future.

###### Acknowledgements.

This work has been supported by the Grant Agency of the Czech Republic (Postdoc project 202/99/D066). We have appreciated critical reading of the manuscript by J. Petzelt from Institute of Physics AS CR.


format: text

text:
# Untitled Document

Origin of Radially Increasing Stellar Scaleheight

in a Galactic Disk

Chaitra A. Narayan and Chanda J. Jog

Department of Physics, Indian Institute of Science

Bangalore 560 012, INDIA.

email: chaitra@physics.iisc.ernet.in,cjjog@physics.iisc.ernet.in

Astron. & Astrophys. Letters, 2002, in press.

Abstract

For the past twenty years, it has been accepted that the vertical scaleheight of the stellar disk in spiral galaxies is constant with radius. However, there is no clear physical explanation for this in the literature. Here we calculate the vertical stellar scaleheight for a self-gravitating stellar disk including the additional gravitational force of the HI and H2 gas and the dark matter halo. We apply our model to two edge-on galaxies, NGC 891 and NGC 4565, and find that the resulting scaleheight shows a linear increase of nearly a factor of two within the optical disk for both these galaxies. Interestingly, we show that the observed data when looked at closely, do not imply a constant scaleheight but actually support this moderate flaring in scaleheight.

Running Title : Radially Increasing Stellar Scaleheight

Key Words: galaxies: fundamental parameters - galaxies: individual - NGC 891, NGC 4565 - galaxies: kinematics and dynamics - galaxies: photometry - galaxies: spiral - galaxies: structure

1. Introdution

In the pioneering studies on vertical luminosity distribution of edge-on stellar disks, van der Kruit and Searle (1981a,b) found that the vertical scaleheight of stars is independent of radius in NGC 4244, NGC 5907, NGC 4565 and NGC 891. This was later confirmed for NGC 4244 by Fry et.al. (1999) and for NGC 891 by Kylafis & Bahcall (1987) and also for other galaxies. In fact, the constant stellar scaleheight has become a well-established paradigm. Recent studies, however, have shown (Kent, Dame & Fazio 1991 and de Grijs & Peletier 1997) a moderate increase with radius in stellar scaleheight.

The constancy of scaleheight comes about when the scalelength with which the vertical velocity dispersion falls off exponentially along the radius, $`h_{\mathrm{vel}}`$, is twice the radial density scalelength, $`h_\mathrm{R}`$, in a galactic disk (van der Kruit and Searle 1981a). However, this relation has no theoretical basis and is used only to obtain the flat scaleheight. We note that if the velocity falls off slower than the above rate ($`h_{\mathrm{vel}}>2h_\mathrm{R}`$), then it results in a flaring disk whereas the converse ($`h_{\mathrm{vel}}<2h_\mathrm{R}`$) gives rise to a disk which tapers down in thickness. These results are valid when the density distribution is defined only by the gravitational force of the stellar disk.

In reality, the molecular and atomic hydrogen gas and the dark matter halo exert a considerable gravitational force on the stellar disk and hence should be included in defining the stellar scaleheight. The effect of these components gains importance at large galactocentric radii, where the surface densities of all the components are comparable to each other, whereas in the inner galactic disk the stellar component dominates over the others. We show that it is because of this additional gravitational force, that the stellar scaleheight is not exponentially flaring at large radii but instead it shows a controlled linear variation. Earlier we have proposed (Narayan & Jog 2002) a gravitationally coupled, three-component disk plus halo model which predicts the scaleheight curves of HI, H2 and the stellar disk of our Galaxy, which are shown to be in good agreement with observations.

In this Letter, we apply the same model for two external galaxies, NGC 891 and NGC 4565, and compare the stellar scaleheights obtained with observations available in the literature. These galaxies are chosen because these are two of the prototypical edge-on galaxies from the sample of van der Kruit & Searle (1981a,b) which also have detailed HI and H2 surface density values from observations.

In Sect. 2, we discuss the details of the model used. The resulting scaleheights for NGC 891 and NGC 4565 are compared with observations in Sect. 3. Section 4 contains discussion and Sect. 5 gives a summary of the conclusions.

2. Calculation for vertical stellar scaleheight

The treatment for obtaining the vertical stellar scaleheight is based on a three-component disk plus halo model, as described by Narayan & Jog (2002), which is briefly summarized below.

The force equation along the $`z`$-axis or the equation of hydrostatic equilibrium is given by (e.g., Rolhfs 1977):

$$\frac{(v_\mathrm{z})_\mathrm{i}^2}{\rho _\mathrm{i}}\frac{\mathrm{d}\rho _\mathrm{i}}{\mathrm{d}z}=(K_\mathrm{z})_\mathrm{s}+(K_\mathrm{z})_{\mathrm{HI}}+(K_\mathrm{z})_{\mathrm{H}_2}+(K_\mathrm{z})_{\mathrm{DM}}$$
$`(1)`$
where $`\rho `$ is the mass density, $`(K_\mathrm{z})=\psi /z`$ is the force per unit mass along the $`z`$-axis, $`\psi `$ is the corresponding potential, and the subscript $`i`$ = s, HI, H2 and DM denotes these quantities for stars, HI, H2 and the dark matter halo respectively. We take the root mean square of the vertical velocities of a component $`(v_\mathrm{z})_\mathrm{i}^2^{\frac{1}{2}}`$ or the random velocity dispersion at a radius $`R`$ and treat the component as being isothermal along $`z`$. The right hand side of equation (1) gives the total vertical force due to all the components.

For a thin disk, the joint Poisson equation reduces to :

$$\frac{\mathrm{d}^2\psi _\mathrm{s}}{\mathrm{d}z^2}+\frac{\mathrm{d}^2\psi _{\mathrm{HI}}}{\mathrm{d}z^2}+\frac{\mathrm{d}^2\psi _{\mathrm{H}_2}}{\mathrm{d}z^2}=\mathrm{\hspace{0.25em}4}\pi G\left(\rho _\mathrm{s}+\rho _{\mathrm{HI}}+\rho _{\mathrm{H}_2}\right)$$
$`(2)`$
The coupled equations (1) and (2) are solved for the density $`\rho `$ versus $`z`$, as a boundary value problem using standard numerical techniques (Press et.al. 1986). To obtain the vertical density distribution of a component, we need to know the surface density of the component which is used as a boundary condition, and its vertical velocity dispersion. This results in a sech2-like density profile, and we use its HWHM to define the scaleheight. We have used a single stellar component for simplicity, though recent observations show (de Grijs, Peletier & van der Kruit 1997) a more sharply peaked luminosity profile which may be explained as a sum of several disk components.

3. Results

3.1. NGC 891

3.1.1. Scaleheights from our model

NGC 891 is an edge-on spiral ($`i=89^\mathrm{o}`$) of Hubble type Sb. It has a very large stellar disk with $`h_\mathrm{R}=4.9`$ kpc (van der Kruit & Searle 1981b). We use a central surface density of 700 M pc-2 which is the value given for our Galaxy (Binney & Tremaine 1987), for the sake of simplicity and also because the results from this section do not depend sensitively on the actual choice of the value used.

Bottema (1993) has measured the vertical velocity dispersions for a set of 12 spirals and finds that the central velocity dispersion is about 60 kms-1 for an Sb type galaxy, which is what we use for NGC 891 and NGC 4565. The dispersion falls off exponentially and is generally modeled (Lewis & Freeman 1989, Bottema 1993) with a scalelength, $`h_{\mathrm{vel}}=2h_\mathrm{R}`$ in accordance with the van der Kruit & Searle (1981a) model - see Sect. 1. A unique value of $`h_{\mathrm{vel}}`$ cannot be derived since there are only three data points for NGC 891. Instead, in our work the ratio $`h_{\mathrm{vel}}/h_\mathrm{R}`$ is allowed to take a range of values between 1-4.

The observed surface densities of HI and H2 are taken from Rupen (1991) and Sofue & Nakai (1993) respectively. The gas velocity dispersions in the Galaxy are 8 kms-1 (Spitzer 1978) and 5 kms-1 (Stark 1984 and Clemens 1985), for HI and H2 respectively. These are taken to be the same for all the spiral galaxies, as observed for HI (Lewis 1984) and H2 (Wilson & Scoville 1990) and as argued theoretically for H2 (Jog & Ostriker 1988). Finally, we calculate the gravitational force of the dark matter halo from a standard model of Mera, Chabrier & Schaeffer (1998) with a total halo mass of $``$ 1012 M. These input parameters are used to solve the coupled equations (1) and (2).

Figure 1 shows the scaleheights obtained by the method discussed in Sect. 2 for the ratio $`h_{\mathrm{vel}}/h_\mathrm{R}`$ = 1, 2, 2.5, 3 and 4. This figure shows that in general, the scaleheight is not flat as deduced by van der Kruit & Searle (1981b). The ratio of less than 2 can be ruled out since it gives unphysically tapering scaleheight.

3.1.2. Scaleheights deduced from observed data

Van der Kruit and Searle (1981b) fit their data using an isothermal stellar disk model with the following space-luminosity density profile as applicable for an edge-on system:

$$\mu (R,z)=\mu (0,0)\left(\frac{R}{h_\mathrm{R}}\right)K_1\left(\frac{R}{h_\mathrm{R}}\right)\mathrm{sech}^2\left(\frac{z}{z_\mathrm{o}}\right)$$
$`(3)`$
where $`\mu (0,0)`$ is the space-luminosity density at the centre of the galaxy, $`z_\mathrm{o}`$ represents the scaleheight and $`K_1`$ is the modified Bessel function of the second kind. We note that, by definition, the scaleheight we obtain (see Sect. 2) is not the same as $`z_\mathrm{o}`$ but the two are equivalent. This model gives a single and well-defined sech2 curve on a composite $`z`$-profile, when the $`z_\mathrm{o}`$ is kept constant. Figure 2 shows the composite $`z`$-profile for NGC 891. The solid line results when $`z_\mathrm{o}=1`$ kpc is kept constant with radius. However, we note that the observed data are spread over an interval about this single curve. We show that a similar spread in the calculated values of surface brightness is obtained when $`z_\mathrm{o}`$ is allowed to increase linearly from 0.75 kpc to 1.25 kpc over 4 disk scalelengths or the entire optical disk. This range of $`z_\mathrm{o}`$ is found by trial and error such that the model values match well with observations. This shows that the spread in the observed data allows for a linear increase of a factor of 1.7 in $`𝐳_\mathrm{o}`$. This increase allows us to constrain the value of $`h_{\mathrm{vel}}/h_\mathrm{R}`$ to be between 2-2.5 (see Fig. 1).

3.2. NGC 4565

NGC 4565 is an Sb type spiral galaxy with an inclination of 86o. Figure 3 shows the stellar scaleheight obtained for NGC 4565 using the method described in Sect. 2. The input parameters used here are as follows. The radial scalelength, $`h_\mathrm{R}`$ is 5.5 kpc (van der Kruit & Searle 1981a). The central surface density of stars is taken to be 700 M pc-2 and the central vertical velocity dispersion = 60 kms-1. The choice of above two numbers is justified in Sect. 3.1.1. Surface densities of HI and H2 are taken from Rupen (1991) and Sofue & Nakai (1994) respectively. Since the observations for the stellar velocity dispersion are not available for NGC 4565, we vary the ratio $`h_{\mathrm{vel}}/h_\mathrm{R}`$ within the range of 1-4.

Using the composite $`z`$-profiles (see Sect. 3.1.2), the stellar scaleheight for NGC 4565 was shown to be constant at $`z_\mathrm{o}=0.79`$ kpc (van der Kruit & Searle 1981a). This was later confirmed by Naslund & Jorsater (1997). As in the case of NGC 891, the observed surface brightness values do not give rise to a sharp and single sech2 curve. Figure 4 shows the numerically obtained composite $`z`$-profile, and also the observed data. The noted spread in observed data can be reproduced by increasing the $`z_\mathrm{o}`$ linearly from 0.4 kpc to 1.0 kpc within the optical disk ($``$ 22 kpc). This shows that the observed data of NGC 4565 can hide as much as a factor of 2.5 variation in scaleheight. This allows us to constrain $`h_{\mathrm{vel}}/h_\mathrm{R}`$ to be between 2.5-3 (see Fig. 3).

4. Discussion

1. The results in Sect. 3 show that the observational data can accommodate a substantial radial variation in scaleheight. However, we would like to caution that the observed spread in data need not be entirely due to variation in scaleheight because, observational uncertainties such as the inclination, dust extinction and large measurement errors in photographic plates could also contribute to it. Morrison et al. (1997) conclude from a more complete analysis of the 2-D data that the scaleheight cannot be assigned a unique value, instead it lies within a range. This supports our model.

2. From a recent study of 48 edge-on spirals, de Grijs & Peletier (1997) find that the galaxies show a small radial increase in scaleheight in agreement with our model. They also find that the increase in scaleheight is smaller for later-type galaxies, and have explained this by using the presence of a thick disk. We note that, alternatively, this trend can be explained by our model. This is because there is an increase in the gas mass fraction (Binney & Merrifield 1998) and the dark matter halo fraction (Broeils 1992) for later Hubble types. Therefore, their confining effect would be higher and hence the scaleheight increase would be smaller.

5. Conclusions

We examine the well-established constancy of vertical stellar scaleheight in spiral galaxies, and show that this result is not justified physically nor is it supported by observations. We show for NGC 891 and NGC 4565 that the observational uncertainties can hide as much as a factor of 2 increase in the scaleheight within the optical disk. We show that this variation can be explained theoretically by our model which takes account of the gravitational coupling between stars, gas, and the dark matter halo. Such moderate increase in scaleheight is likely to be common in other spiral galaxies as well.

Acknowledgements We would like to thank the referee, R. de Grijs, for his detailed and constructive comments which have greatly improved the paper.

References

Binney, J., & Tremaine, S. 1987, Galactic Dynamics

(Princeton: Princeton Univ. Press)

Binney, J. & Merrifield, M. 1998, Galactic Astronomy

(Princeton: Princeton Univ. Press)

Bottema, R. 1993, A&A, 275, 16.

Broeils, A.H. 1992, Ph.D. thesis, Groningen Univ.

Clemens, D.P. 1985, ApJ, 295, 422.

de Grijs, R. & Peletier, R.F. 1997, A&A, 320, L21.

de Grijs, R., Peletier, R.F., & van der Kruit, P.C. 1997,

A & A, 327, 966

Fry A.M., Morrison, H.L., Harding, P., & Boroson, T.A.

1999, AJ, 118, 1209.

Jog, C.J. & Ostriker, J.P. 1988, ApJ, 328, 404.

Kent, S.M., Dame, T. & Fazio, G. 1991, ApJ, 378, 131.

Kylafis, N.D., & Bahcall, J.N. 1987, ApJ, 317, 637.

Lewis, B.M. 1984 ApJ, 285, L453.

Lewis, J.R. & Freeman, K.C. 1989, AJ, 97, 139.

Mera, D., Chabrier, G., & Schaeffer, R. 1998, A&A, 330,

953

Morrison, H.L., Miller, E.D., Harding, P., Stinebring, D.R.

& Boroson, T.A. 1997, AJ, 113, 206.

Narayan, C.A. & Jog, C.J. 2002, Submitted to A&A.

Naslund, M. & Jorsater, S. 1997, A&A, 325, 915.

Press, W.H., Flannery, B.P., Teukolsky, S.A., &

Vetterling, W.T. 1986, Numerical Recipes (Cambridge:

Cambridge Univ. Press), chap. 6.

Rohlfs, K. 1977, Lectures on Density Wave Theory, (Berlin

: Springer-Verlag).

Rupen, M.P. 1991, AJ, 103, 48.

Sofue, Y. & Nakai, N. 1993, PASJ 45, 139.

Sofue, Y. & Nakai, N. 1994, PASJ 46, 147.

Spitzer, L. 1978, Physical Processes in the Interstellar

Medium (New York: John Wiley)

Stark, A.A. 1984, ApJ, 281, 624

van der Kruit, P.C. & Searle, L. 1981a, A&A, 95, 105.

van der Kruit, P.C. & Searle L. 1981b, A&A, 95, 116.

Wilson, C.D. & Scoville, 1990, ApJ, 363, 435.

Figure Legends

Figure 1. Variation in resulting stellar vertical scaleheight with galactocentric radius for NGC 891. The variation is found to crucially depend on the ratio $`h_{\mathrm{vel}}/h_\mathrm{R}`$. While the scaleheight is found to decrease by a factor of 50 and 1.1 for a ratio of 1 and 2 respectively, and increase by a factor of 1.8, 2.8 and 4.5 for the ratio of 2.5, 3 and 4 respectively, within the optical disk ($``$ 20 kpc). Thus the scaleheight need not be constant with radius for the galaxy.

Figure 2. The plot of surface brightness, $`\mu `$, versus the distance from the midplane, $`z`$, for NGC 891. This composite $`z`$-profile is obtained by vertically shifting the individual $`z`$-profiles into coincidence at $`z=1.5`$ kpc. The solid line results when $`z_\mathrm{o}=1`$ kpc and is independent of radius. When $`z_\mathrm{o}`$ is linearly increased from 0.75-1.25 kpc within 20 kpc, surface brightness values (points) disperse over an interval of 1 mag/arcsec2 at $`z=0`$, as is exactly seen in the original data (van der Kruit & Searle 1981b). The vertical bars indicate the range of the data. This shows that the data can allow for as much as a factor of 1.7 increase in scaleheight. The deviation of the observed data at high $`z`$ is due to the thick disk which is not included in the model.

Figure 3. Variation in resulting stellar vertical scaleheight with radius for NGC 4565. The scaleheight is found to decrease by a factor of 50 and 1.1 for a ratio of 1 and 2 respectively, and increase by a factor of 1.8, 2.7 and 4.4 for the ratio of 2.5, 3 and 4 respectively, within the optical disk ($``$ 22 kpc). Thus the scaleheight need not be constant with radius for this galaxy too.

Figure 4. The plot of surface brightness, $`\mu `$, versus the distance from the midplane, $`z`$, for NGC 4565, with details similar to Figure 2. The linear increase in $`z_\mathrm{o}`$ from 0.4-1.0 kpc or an increase of a factor of 2.5 in the scaleheight can be easily hidden by the spread in the observed data.


text:
# A slow and dark atomic beam

## 1 Introduction

Cold atomic beams are useful for atom interferometry and precision atomic spectroscopy experiments. Cold beams were first derived from oven sources and slowed by the radiation pressure from a resonant laser beam which can be temporally chirped  or operated at fixed frequency in conjunction with Zeeman coils . These atomic beams have a large transverse velocity spread because the atoms are cooled in only one dimension. Colder atomic beams that are both transversely and longitudinally cooled can be derived from 2D magneto-optical traps (MOTs) and moving optical molasses . The mean beam velocity produced by these 2D MOTs is tunable between 0.5-3 m/s. Atoms can also be extracted from 3D MOTs using the Low Velocity Intense Source (LVIS) geometry, in which a radiation pressure imbalance is created by a small hole at the center of one of the retro-reflecting mirror-waveplate assemblies , or by an additional laser beam . These sources produce atomic beams with a mean velocity of about 15 m/s. Recently, the propagation of cold atomic beams in magnetic guides has been studied .

Evaporative cooling of transversely confined atomic beams is a promising route towards the generation of a continuous source of Bose Einstein Condensates (BEC) . This approach to a continuous coherent atom source differs from the method of periodic replenishment of a static BEC reservoir , where suitably interlaced moving optical dipole traps are used. In the latter method, large number and phase fluctuation of the condensate can be expected. The phase and atom number will be more stable if the whole system runs continuously, without moving potentials and/or pulsed atomic beams. A BEC continuously replenished by an evaporatively cooled, guided atomic beam should therefore have much lower phase and atom number variations.

Evaporative cooling of a magnetically guided atomic beam presents technical challenges. To achieve high enough densities and collision rates, a very large flux of cold atoms (better than $`10^9`$ s-1) moving at very low velocities (less than 1 m/s) is required. The atomic flux must be in a low-field-seeking state and must be efficiently coupled into the guide. These challenges can only be partially met when using the beam sources described in , which produce atomic beams with longitudinal velocities of order 15 m/s. Moving optical molasses combined with magneto-optic trapping can produce dense, guided atomic beams with velocities in the 1 m/s range . Other schemes to couple cold atoms into magnetic guides have been studied in .

The simultaneous operation of a MOT producing an atomic beam and a nearby magnetic atomic-beam guide is a further challenge. In such a system, the MOT will produce some stray light in its vicinity. The MOT stray light will be scattered by the atoms propagating in the atom guide, resulting in optical pumping of the atoms out of the desired low-field-seeking magnetic sublevel. Stray light intensities $`<1`$ nW/cm2 of near-resonant light (for Rb atoms) will seriously attenuate the magnetically guided atomic beam. It appears very hard to suppress the MOT stray light below that intensity level. A promising way out of this dilemma is to optically pump the atomic beam emitted from the MOT into a dark, low-field-seeking state before the near-resonant stray light from the MOT depletes it. In alkali atoms, a magnetic sublevel of the lower of the two ground-state hyperfine levels can be used as such a dark state. To ensure a rapid enough transfer into that state, a separate optical-pumping laser can be used. In this paper, we use this idea to produce a slow beam of 87Rb atoms in the $`|F=1,m_\mathrm{F}=1`$ ground-state sublevel, which is dark with respect to the stray light of a nearby MOT (which is the atomic-beam source). In addition to being in a dark state, the atomic beam is also considerably slower than the bright atomic beams in . Experimental data from pulsed measurements are shown, but the method can also generate continuous, dark atomic beams. The beam may be magnetically guided and, in future efforts, evaporatively cooled while the MOT is operating.

## 2 Extraction of Atoms from a Vapour Cell MOT

We use a six-beam vapour-cell MOT and a push laser to produce a 87Rb atomic beam. The radiation pressure imbalance due to the push beam  accelerates cold atoms out of the center of the MOT. The MOT and pusher beams act on the upper hyperfine level $`F=2`$ of 87Rb. The velocity of the resultant atomic beam is determined by the intensity, detuning, and polarization of the push laser. Since the Zeeman and the Doppler shifts of the atoms are position- and/or velocity-dependent, the acceleration also depends on position and velocity. On average, a 87Rb atom extracted from the MOT scatters about 10,000 photons from the push beam before it is pumped into the lower hyperfine state $`F=1`$ via an accidental off-resonant $`F=2F^{}=2`$ transition. The acceleration phase of the atom can, however, be terminated early using an additional laser beam, referred to as transfer beam, that optically pumps the atom into the $`F=1`$ ground state. Because of the large energy spacing between the $`F=1`$ and $`F=2`$ ground states, the push and the MOT light is not resonant with the $`F=1`$ atoms. The dark atomic beam is slower than that produced by a conventional LVIS. Also, all atoms in the beam cease accelerating at the same location, namely when they enter the transfer beam. This fact leads to a velocity spread of the dark atomic beam that is less than the velocity spread of a conventional LVIS atomic beam.

The repumper light required to operate the MOT is resonant with the $`F=1`$ $``$ $`F^{}=2`$ transition and must therefore be sufficiently shielded from the $`F=1`$ dark atomic beam. There is only a single MOT repumper beam, and that beam has a low intensity and is quite small. Therefore, it is technically much easier to shield the repumper stray light from the dark atomic beam than preventing the MOT stray light, which is near-resonant with the $`F=2`$ $``$ $`F^{}=3`$ transition, from reaching the atomic beam. This (practical) difference in the ability to shield MOT stray light from the atomic beam vs. the ability to shield repumper stray light has been a major motivation for the development of the method presented in this paper.

## 3 Experimental Setup

Our experimental setup has been described previously in  and is shown schematically in Fig. 1. We use a two-wire magnetic atom guide. The guide wires are separated by 2.8 cm and are placed outside a quartz vacuum tube. The transverse magnetic field of the atom guide varies along the guide axis ($`z`$ direction) due to the presence of two tapered pieces of magnetic iron. At a current of 400 A, the transverse magnetic field gradient increases from 70 G/cm to 320 G/cm along the guide. 87Rb atoms are captured from the background gas in a MOT, which is formed on the axis of the guide. The MOT laser beams are detuned by a frequency of $`2\mathrm{\Gamma }/2\pi =12`$ MHz to the red of the $`5S_{1/2}F=25P_{3/2}F^{}=3`$ transition and have an intensity of 5 mW/cm2. An additional beam resonant with the $`5S_{1/2}F=15P_{3/2}F^{}=2`$ repumps atoms that fall into the $`F=1`$ state back into $`F=2`$. In our pulsed experiments, during the MOT loading phase the guide current is reduced to 200 A, and a pair of external magnetic coils is energized, resulting in a MOT magnetic field with gradients -15, 26 and -11 G/cm in the $`x,y`$ and $`z`$ directions. In each loading cycle, $`2\times 10^7`$ atoms are collected by the MOT.

Cold atoms are extracted from the MOT using a push laser beam propagating along the guide axis . In the presented experiments, the guide current (see Fig. 1) is increased from 200 A to 400 A at the time when the push laser beam turns on, while the current in the external coils is kept constant. The push beam has the same detuning as the six MOT beams and produces a radiation pressure imbalance at the MOT center that forces the atoms out. The MOT repumper beam intersects the $`z`$-axis at a right angle and is masked by a sharp edge such that the repumper intensity abruptly drops to zero a few mm downstream from the MOT (see Fig. 1). The transfer beam, located close to the MOT, optically pumps the extracted atoms into the $`F=1`$ dark state. The transfer beam is locked to the $`5S_{1/2}F=25P_{3/2}F^{}=2`$ transition. The transfer beam has a circular cross section (diameter of $`2`$ mm) with an approximately constant intensity equal to the saturation intensity ($`I_{\mathrm{sat}}=1.6`$ mW/cm2). The transfer beam interrupts the acceleration of the atoms. While providing a well defined location of transfer from $`F=2`$ to $`F=1`$, repeated transitions of the atoms back and forth between these states must be avoided. Therefore, any spatial overlap between the repumper and the transfer beam needs to be avoided.

In our experiments we have, for convenience, used an existing small uncoated quartz cell as our vacuum chamber. Therefore, the amount of stray light present in the vacuum system has been unusually high. In order to reduce the amount of repumper stray light in the atom guiding region to an acceptable level, the MOT repumper intensity had to be reduced to levels so low that the number of atoms in the MOT has been significantly reduced. Obviously, a larger vacuum chamber with dark surfaces and vacuum windows placed outside direct line-of-sight from the guided atoms will allow us to use a higher MOT repumper intensity.

## 4 Experimental Data

To detect the fluorescence of the atomic beam, the beam is imaged onto a 3 mm by 3 mm photodiode (magnification approximately one). The center of the photodiode corresponds to an object location that has a distance of 1 cm from the MOT center. To measure the time-of-flight distribution of the atomic beam, a mechanical shutter is used to pulse the push beam. The onset of the push beam extracts a pulse of atoms from the MOT. The MOT beams can be turned off or left on continuously during the extraction. The extracted atoms stream through the region that is imaged onto the photodiode. Atoms that enter the detection region in the state $`F=2`$ are still being accelerated by the push beam. The resultant fluorescence is detected, and its time-dependence can be used as a measure for the average velocity and velocity spread of the atoms. To be able to also detect atoms propagating in the dark state ($`F=1`$), an additional repumping beam can be introduced into the detection region. This repumper returns atoms propagating in the dark state $`F=1`$ into the bright state $`F=2`$, thereby making them detectable.

### 4.1 Pulsed Measurements of Velocity

The following data shows the time dependence of the fluorescence signal for the case that the leading edge of the transfer beam is located 3 mm from the MOT. The MOT beams are shut off before the push beam is turned on. Without the transfer beam, most of the atoms remain in the $`F=2`$ bright state and scatter light from the push beam (solid curve in Fig. 2) while accelerating. Assuming constant acceleration we find, based on the average arrival time of the atoms, an average velocity of the atoms of about 10 m/s. When the transfer beam is turned on and the additional repumper in the detection region is left off, most of the atoms are pumped into the $`F=1`$ dark state and traverse the detection region undetected. The weak remnant signal in Fig. 2 corresponds to a small fraction of atoms that are still in the bright state $`F=2`$. The presence of those atoms is due to stray repumper light throughout the chamber that accidentally returns some atoms from the dark into the bright state before they leave the detection region.

The dashed curve in Fig. 2 shows the signal observed when both the transfer beam and the additional repumper beam in the detection region are used. In that case, the acceleration of the atoms is suspended at the leading edge of the transfer beam, and resumes at the leading edge of the repumper beam in the detection region. The fluorescence associated with the renewed acceleration is measured. The peak of this signal is delayed by 0.4 ms with respect to the signal measured without transfer beam. Assuming that the transfer beam is located approximately 3 mm from the MOT, from the observed delay of 0.4 ms we deduce an average velocity of the dark beam of about 5 m/s. This is more than a factor of two slower than any atomic beam based on the LVIS design. The velocity of the dark atomic beam can be controlled by the intensity of the push beam as well as the distance of the transfer beam from the MOT center. We demonstrate this in the next two sections.

### 4.2 Effect of the Push Beam Intensity

In the following, we discuss a series of fluorescence measurements obtained for various push beam intensities. The rising edge of the transfer beam is located 3 mm from the MOT center. The MOT beams are always on, and the radiative force of the pulsed push beam must overcome the damping force of the MOT to extract the atoms. We find that for a given push beam intensity the push beam is most effective if it is $`\sigma `$-polarized, with a helicity that is opposite to that one would use if the pusher beam was a MOT trap beam. Under this condition, the magnetic field of the MOT tunes the atoms closer to resonance with the push beam as they are pushed out from the MOT, resulting in a larger acceleration. We have also found that for the push beam to overcome the damping force of the MOT, the push beam power must exceed a threshold power of about 0.1 mW. Below this power, the push beam merely displaces the MOT, and there is no atomic beam.

Fig. 3(a) shows the fluorescence signal for different values of the push beam power when both the transfer beam and the additional repumper laser in the detection region are turned on. As the push beam power decreases, the fluorescence signal peaks at later times, corresponding to a lower average beam velocity. For each push beam power, we have, for comparison, also measured the fluorescence signals without transfer laser (data not shown). For the resultant pairs of curves, one pair being analogous to the solid and the dashed curve of Fig. 1, we have then determined the difference in the most probable arrival time of the atoms. In Fig. 3(b), the time differences are plotted vs. push beam power. With decreasing power, the delay time between the transfer-on and transfer-off signals increases from 0.02 ms to 1.00 ms. This behavior reflects the fact that the time of travel of the dark atomic beam from the transfer beam to the detection region increases with decreasing velocity.

### 4.3 Effect of the Transfer Beam Position

The velocity of the dark atomic beam can also be controlled by varying the position of the transfer beam relative to the MOT center. Moving the transfer laser closer to the MOT leaves the atoms less time to interact with the push beam before they are pumped into the dark state. This results in a slower dark atomic beam. Fig. 4 shows the fluorescence signal, with both the transfer beam and the additional repumper beam on, for three distances of the transfer beam from the MOT. As explained in Sec. 3 and in Fig. 1, the cutoff of the MOT repumper beam and the transfer beam need to be varied such that there always remains a small gap between both beams. As expected, the delay time caused by the transfer beam increases when the transfer beam is moved towards the MOT. There is a tradeoff between attaining a large atom flux and a low atomic-beam velocity. While moving the transfer laser beam closer to the MOT results in a slower atomic beam (which may be desired), it also entails the necessity to move the MOT repumper cutoff closer to the MOT center, thereby decreasing the MOT loading efficiency and the atomic-beam flux.

## 5 Numerical Simulation

We have numerically simulated the dynamics of the MOT with the push and transfer beams using a classical description of the center-of-mass motion and a quantum-mechanical description of the internal dynamics of the atoms. We have used the analysis given in  as a starting point for our model. For a two-level system, the radiation pressure force acting on the atom due to a single laser beam is the product of the photon momentum multiplied by the photon scattering rate $`\mathrm{\Gamma }_{\mathrm{scatter}}`$, and is given by

$`𝐅=\mathrm{}𝐤\mathrm{\Gamma }_{\mathrm{scatter}}=\mathrm{}𝐤{\displaystyle \frac{\mathrm{\Gamma }}{2}}{\displaystyle \frac{s}{s+1}},`$ (1)
where $`\mathrm{\Gamma }`$ is the decay rate of the excited state and $`s`$ is the saturation parameter

$`s={\displaystyle \frac{I}{I_{\mathrm{sat}}}}{\displaystyle \frac{\mathrm{\Gamma }^2}{\mathrm{\Gamma }^2+4(\delta 𝐤𝐯\mathrm{\Delta }E_{\mathrm{Zeeman}}/\mathrm{})^2}}.`$ (2)
The saturation parameter takes into account the Doppler shift due to the (classical) atomic velocity $`𝐯`$ as well as the relative Zeeman shift $`\mathrm{\Delta }E_{\mathrm{Zeeman}}`$ between the involved states arising from the magnetic field of the MOT coils and the atom guide wires. Further, $`\delta `$ is the field-free laser detuning, $`I`$ the beam intensity, and $`I_{\mathrm{sat}}=1.6`$ mW/cm2 the saturation intensity.

The $`F=2F^{}=3`$ MOT transition of 87Rb involves multiple magnetic subcomponents with respective Zeeman shifts and Clebsch-Gordan coefficients. Further, to be able to calculate the average radiation pressure force, the probability distribution of the atom over the magnetic sublevels of the $`F=2`$ ground state needs to be determined. To make the model numerically more tractable, we use an $`F=0F^{}=1`$ model transition with wavelength and saturation intensity equal to those of the Rb D2 line, as explained in . Also, due to the low density of atoms in MOTs designed to produce cold atomic beams, we can ignore collisions between atoms and simulate the dynamics of each particle independently.

Using these assumptions, the radiation-pressure force on an atom due to the MOT beams and the pusher beam becomes

$`𝐅={\displaystyle \underset{i=1,7}{}}{\displaystyle \underset{m=1,1}{}}\mathrm{}𝐤_\mathrm{i}{\displaystyle \frac{\mathrm{\Gamma }}{2}}{\displaystyle \frac{s_{\mathrm{i},\mathrm{m}}}{1+\underset{i^{},m^{}}{}s_{\mathrm{i}^{},\mathrm{m}^{}}}}`$ (3)
where

$`s_{\mathrm{i},\mathrm{m}}={\displaystyle \frac{I_{\mathrm{i},\mathrm{m}}}{I_{\mathrm{sat}}}}{\displaystyle \frac{\mathrm{\Gamma }^2}{\mathrm{\Gamma }^2+4(\delta _\mathrm{i}𝐤_\mathrm{i}𝐯+m\mu _\mathrm{B}/\mathrm{})^2}}`$ (4)
is the saturation parameter of the individual beam components. The $`i`$ and $`m`$ are indices for the beam number (six MOT beams plus one push beam) and the laser polarization, respectively. Each laser beam $`i`$ is decomposed into its $`\sigma _{}`$, $`\pi `$ and $`\sigma _+`$ intensity components $`I_{\mathrm{i},1},I_{\mathrm{i},0}`$, and $`I_{\mathrm{i},1}`$ in a coordinate frame aligned with the magnetic-field direction at the (time-dependent) location of the atom. The denominator term $`1+\underset{i^{},m^{}}{}s_{\mathrm{i}^{},\mathrm{m}^{}}`$ in Eq. 3 is introduced in order to obtain a reasonable saturation behavior  (the magnitude of the force on an atom must not exceed $`\mathrm{}k\mathrm{\Gamma }/2`$).

In the following, we explain the features we have added to the model outlined in . Spontaneous emission causes a diffusive motion in momentum space. In our model, we integrate, as the atom propagates, the total number of photons $`N`$ that have been scattered. At suitable time intervals, a velocity kick of magnitude $`\sqrt{N}`$ times the single-photon recoil velocity and with a random direction is added to the velocity of the atom, and $`N`$ is re-set to zero. We have verified that this procedure produces the expected steady-state temperature and cloud size of the MOT (with the push beam turned off). Further, in addition to the ground and excited states, we introduce a virtual dark state that the atom can be transferred to by a transfer beam. This state does not interact with the MOT or push beams. The role of the repumper is simulated by a transfer of atoms back to the $`F=0`$ state of the simulation. We ignore the momentum changes associated with the repumper and transfer beams because they amount to just a few photon recoils.

Fig. 5(a) shows the simulated fluorescence signals corresponding to the conditions for the data shown in Fig. 4. The calculations reproduce the qualitative features of the effect of the push beam position. The corresponding velocity distributions of the atoms are shown in Fig. 5(b). The simulations indicate a much lower beam velocity than the experimental data. This quantitative discrepancy can be attributed in part to imperfect beam geometries in the experiment, and, mostly, to repumper stray light scattered by the uncoated surfaces of our quartz vacuum chamber. The repumper stray light can re-pump atoms that have been transferred into the dark state back into the bright state and thereby make them susceptible to unwanted acceleration.

The atomic beam velocity is sensitive to all laser beam parameters. In our simulation, suitable parameters lead to dark atomic beams with velocities as slow as 0.3 m/s, which is about a factor 30 smaller than the method of the LVIS. A few simulated velocity distributions are shown in Fig. 6. There, the RMS velocity spread is of order 0.2 m/s in all cases. The slowest velocites are obtained using push beam intensity of 3 mW/cm2 and a transfer beam separation from the MOT of 3 mm. While experimental limitations have prevented us from a quantitative reproduction of these theoretical results, we have been able to demonstrate the basic functionality of the technique.

## 6 Conclusion

We have experimentally demonstrated a simple technique to produce atomic beams with velocities and velocity spreads that are much smaller than that obtained from an LVIS. By pumping the atoms into a dark state, we have produced a slow, dark atomic beam that can be guided in a magnetic atom guide. Since the guided, dark atomic beam is not susceptible to stray light from the nearby MOT atom source, it can be operated continuously and concurrently with the MOT. Our numerical simulations show that the technique can produce dark atomic beams with velocities less than 0.2 m/s and velocity spreads comparable to that of a (doppler-limited) moving optical molasses. Our small vacuum chamber with uncoated windows has resulted in significant stray repumper light in the guide, which prevented us from observing the very slow velocities predicted by our simulations. This problem can, however, easily be alleviated using a larger vacuum chamber with fewer reflecting surfaces. We believe that dark, guided atomic beams are an ideal atom source for the realization of a continuous BEC in a magnetic atom guide.


text:
# Parallel versus sequential updating for Belief Propagation decoding

## I Introduction

Error correcting codes are essential part of modern communication, enabling reliable transmission over noisy channels. Bounds over the channels capacity were derived by Shannon in 1948 Shannon , but more than four decades passed before codes that nearly saturate the bound, such as Turbo code Turbo  and LDPC shokrolahi spilman irreg 65-85  were presented. In recent years, an interesting bridge was established between error correcting codes and statistical mechanics of disordered systems Sourlas ; sourlas - montanari ; saad -kabashiba: SM of ECC .

It is well known that as the noise level in a channel increases, the decoding time (measured in algorithm iterations) also increases.(shokrollahi-richrdson-urbanke 3, , section 3.2). Furthermore, as the noise $`f`$ approaches the threshold level, $`f_c`$, the number of iterations diverges as a power-law, $`t1/(f_cf)`$kanter saad construction . Hence, the acceleration of the decoding process, or the reduction of the required number of iterations, is essential to ensure a smooth information flow when operating near the channels capacity.

In this paper we propose a variation of the well known Belief Propagation (BP) algorithm mackay-neal-short-article4 ; Gallager ; Mackey good error 1 , which we label as "Sequential Updating Scheme" (SUS). The complexity per iteration of the SUS is similar to the complexity per iteration of the *regular* Belief Propagation (BP) algorithm, with Parallel Updating Scheme (PUS). However, simulations over a Binary Symmetric Channel (BSC) indicate that for a given code, the SUS requires about *one half* of the iterations in comparison to the PUS, while the averaged bit error probability, $`p_b`$, is the same.

This article is organized as follows: In section 2 the parallel and the sequential updating schemes are defined. The distribution of the decoding time obtained in simulations over BSC for both schemes are presented in section 3. In section 4 the complexities of the sequential and the parallel updating schemes are compared. A qualitative theoretical argument supporting the acceleration of the decoding procedure in the SUS is presented in section 5, and a detailed description of the simulations is presented in section 6. A brief conclusion is presented in section 7.

## II The BP algorithm

### II.1 Notation

We focus on transmission over a BSC where each transmitted bit has a chance $`f`$ to flip during transmission, and a chance $`1f`$ of being transmitted correctly. Most of our results were obtained using Mackey and Neal’s algorithm (MN) , described in detail in reference Mackey good error 1 . Briefly, the algorithm is defined as follows:

A word of size $`n`$ is encoded into a codeword of size $`m`$, ($`Rate=n/m`$), using the following binary matrices:

a sparse matrix of dimensions $`(m\times n)`$

a sparse and invertible matrix of dimensions $`(m\times m)`$.

The encoding of a word $`s`$ into a codeword $`t`$ is performed by:

$$t=B^1As(mod\mathrm{\hspace{0.25em}2})$$
(1)
During the transmission, a noise $`n`$ is added to the data, and the received codeword $`r`$ is:

$$r=t+n(mod\mathrm{\hspace{0.17em}2})$$
(2)
The decoding is performed by calculating $`z=Br`$

$$z=Br=B(t+n)=B(B^1As+n)=As+Bn=[A,B][s,n],$$
(3)
where \[ \] represents appending matrices or concatenating vectors.

Denoting $`H=[A,B]`$ and $`x=[s,n]`$, the decoding problem is to find the most probable $`x`$ satisfying: $`Hx=z`$ (mod2) , where:

is $`(m\times (n+m))`$ matrix

is the constraints (checks) vector of size $`m`$.

is the unknown (variable) vector of size $`n+m`$, termed variable vector.

This problem can be solved by a BP algorithm. Following Mackey good error 1 ; kanter - saad cascading 7 , we refer to the elements of $`x`$ as nodes on a graph represented by $`H`$, the elements of $`z`$ being values of checks along the graph. The non-zero elements in a row $`i`$ of $`H`$ represent the bits of $`x`$ participating in the corresponding check $`z_i`$. The non-zero elements in column $`j`$ represent the checks in which the $`j`$th bit participates.

We follow Kanter and Saad’s (KS) construction for the matrices $`A,B`$ KS gaussian ; kanter - saad cascading 7 ; kanter saad construction . This construction is characterized by very sparse matrices and a cyclic form for $`B`$, together with relatively high error correction performance.

For each non-zero element in $`H`$, the algorithm calculates 4 values mackay-neal-short-article4 . The coefficient $`q_{ij}^1`$ ($`q_{ij}^0`$) stands for the probability that the bit $`x_j`$ is 1 (0), taking into account the information of all checks in which it participates, except for the $`i`$th check. The coefficient $`r_{ij}^1`$ ($`r_{ij}^0`$) indicates the probability that the bit $`x_j`$ is 1 (0), taking into account the information of all bits participating in the $`i`$th check, except for the $`j`$th bit.

The algorithm is initialized as follows. The coefficient $`Q_j`$ is set equal to our prior knowledge about that bit $`j`$. In our simulations we assume $`Q_j=0.5`$ if it is a source bit ($`jN`$), and $`Q_j=f`$ if it is a noise bit ($`j>N`$). Then the $`q`$ values are set: $`q_{ij}^1=Q_j`$ ; $`q_{ij}^0=1q_{ij}^1`$ for all non-zero elements in the $`j`$th column.

### II.2 Parallel Updating Scheme (PUS)

The PUS consists of alternating horizontal and vertical passes over the $`H`$ matrix. Each pair of horizontal and vertical passes is defined as an iteration. In the horizontal pass, all the $`r_{ij}`$ coefficients are updated, row after row:

$$r_{ij}^0=\underset{(allconfigurationswithx_j=0,satisfingz_i)}{}\underset{j^{}j}{}q_{ij^{}}^{x_j^{}}$$
(4)
$$r_{ij}^1=\underset{(allconfigurationswithx_j=1,satisfingz_i)}{}\underset{j^{}j}{}q_{ij^{}}^{x_j^{}}$$
(5)
where it is clear that the multiplication is performed only over the non-zero elements of the matrix $`H`$.

A practical implementation of (4) and (5) is carried out by computing the differences $`\delta q_{ij}q_{ij}^0q_{ij}^1`$, and $`\delta r_{ij}r_{ij}^0r_{ij}^1`$ is then obtained from the identity:

$$\delta r_{ij}=(1)^{z_i}\underset{j^{}j}{}\delta q_{ij^{}}.$$
(6)
From the normalization condition $`r_{ij}^0+r_{ij}^1=1`$ one can find:

$$r_{ij}^0=(1+\delta r_{ij})/2;r_{ij}^1=(1\delta r_{ij})/2$$
(7)
(For a detailed description of this method, see Mackey good error 1 ).

In the vertical pass, all $`q_{ij}^1,q_{ij}^0`$ are computed, column by column, using the updated values of $`r_{ij}^1,r_{ij}^0`$:

$$q_{ij}^0=\alpha _{ij}p_j^0\underset{i^{}i}{}r_{i^{}j}^0$$
(8)
$$q_{ij}^1=\alpha _{ij}p_j^1\underset{i^{}i}{}r_{i^{}j}^1$$
(9)
where $`\alpha _{ij}`$ is a normalization factor, chosen to satisfy $`q_{ij}^0+q_{ij}^1=1`$, and $`p_j^0`$ , $`p_j^1`$ are the priors. Now the pseudo-posterior probability vector $`Q`$ can be computed by:

$$Q_j^0=\alpha _jp_j^0\underset{i}{}r_{ij}^0$$
(10)
$$Q_j^1=\alpha _jp_j^1\underset{i}{}r_{ij}^1$$
(11)
Again, $`\alpha _j`$ is a normalization constant satisfying $`Q_j^1+Q_j^0=1`$ , and $`i`$ runs only over non-zero elements of $`H`$. Each iteration ends with a clipping of $`Q`$ to the variable vector $`x`$ : if $`Q_j>0.5`$ then $`x_j=1`$, else: $`x_j=0`$. At the end of each iteration a convergence test, checking if $`x`$ solves $`Hx=z`$, is performed. If some of the $`m`$ equations are violated, the algorithm turns to the next iteration until a pre-defined maximal number of iterations is reached with no convergence (our halting criteria are described in detail in section 6). Note that there is no inter-iteration information exchange between the bits: all $`r_{ij}`$ values are updated using the previous iteration data.

### II.3 Sequential Updating Scheme (SUS)

In the SUS, we perform the horizontal and vertical passes separately for each bit in $`x`$. A single sequential iteration for the bit $`x_j`$ consists of the following steps:

1. For a given $`j`$ all $`r_{ij}`$ are updated. More precisely, for all non-zero elements in column $`j`$ of $`H`$, use (6,7) for updating $`r_{ij}`$. Note that this is only a *partial* horizontal pass, since only $`r_{ij}`$’s belonging to a specific column are updated.
2. After all $`r_{ij}`$ belonging to a column $`j`$ are updated, a vertical pass as defined in (8,9) is performed over column $`j`$. Again, this is a *partial* vertical pass, referring only to one column.
3. Steps 1-2 are repeated for the next column , until all columns in $`H`$ are updated.
4. Finally, the pseudo posterior probability value $`Q_j`$, is calculated by (10,11).

After all variable nodes are updated, the algorithm continues as for the parallel scheme: clipping, checking the validity of the $`m`$ equations and proceeding to the next iteration.

## III Results

We performed simulations of decoding over a BSC using various rates, block length and flip rates ($`f`$) and with the following constructions: (a) KS kanter saad construction  construction; and (b) Irregular LDPC codes, following the Luby - Mitzenmacher - Shokrollahi - Spielman construction (LMSS), described in shokrolahi spilman irreg 65-85 . Note that the LMSS code differs slightly from the MN code, but still has a similar form of the BP scheme. We compare the distribution of the convergence times of the PUS and SUS by decoding of the same codewords (samples).

Figure 1 presents the distribution of the convergence times (measured in algorithm iterations) for PUS (filled bars) and SUS (empty bars) . The code rate is 1/2, $`f=0.08`$ (the channel capacity is $`0.11`$ ) and the block length is $`10,000`$. The statistics were collected over at least $`3,000`$ different samples. The converging time for the SUS is about one half of the converging time for the PUS. The average convergence time for the PUS is 32.12 iterations for the KS construction; 28.52 for the LMSS construction; while for the SUS the average convergence time is 16.7 and 16.32 iterations, respectively. It is worth mentioning that the superior decoding time does not damage the error correcting property of the code. In both constructions the observed bit error rate, $`p_b`$, after PUS or SUS, is nearly the same (see Table 1 for details).

In Figure 2 the ratio between the converging times, (SUS/PUS) per sample is plotted (time is measured in algorithm iterations). For the vast majority of the samples this is very close to the average rate. This indicates that the double number of iterations for the PUS in comparison to the SUS is the typical result.

Table 1 presents similar measurements for other rates and noise levels. Results indicate the following general rule. Independent of the construction, the noise level and the rate, the convergence time of the PUS is around double the number of iterations required to achieve convergence in the SUS. (As our statistics was collected over $`3000`$ samples of block size $`10^4`$, we do not report the exact value for $`p_b10^5`$ ).

| construction | $`\begin{array}{c}<t_{PUS}>\\ (iterations)\end{array}`$ | $`\begin{array}{c}<t_{SUS}>\\ (iterations)\end{array}`$ | $`\frac{t_{SUS}}{t_{PUS}}`$ | $`\frac{t_{SUS}}{t_{PUS}}`$ | bit error rate (PUS, SUS) |
| --- | --- | --- | --- | --- | --- |
| KS,R=1/2, $`f`$=0.09 | 58.3 | 28.4 | 0.488 | 0.504 | $`5.410^5`$, $`5.510^5`$ |
| KS, R=1/2, $`f`$=0.08 | 32.12 | 16.7 | 0.519 | 0.5076 | $`p_b10^5,p_b10^5`$ |
| KS, R=1/2, $`f`$=0.07 | 23.0 | 11.54 | 0.505 | 0.501 | $`p_b10^5,p^b10^5`$ |
| KS, R=1/3, $`f`$=0.159 | 65.58 | 32.54 | 0.496 | 0.5236 | $`3.410^4,2.0810^4`$ |
| KS,R=1/3, $`f`$=0.15 | 40.53 | 21.28 | 0.525 | 0.525 | $`p_b10^5,p_b10^5`$ |
| KS, R=1/5, $`f`$=0.23 | 103.9 | 54.66 | 0.526 | 0.5005 | $`1.910^3,1.910^3`$ |
| KS, R=1/5, $`f`$=0.21 | 36.98 | 19.46 | 0.527 | 0.526 | $`p_b10^5,p_b10^5`$ |
| LMSS, R=1/2, $`f`$=0.08 | 28.52 | 16.32 | 0.572 | 0.573 | $`p_b10^5,p_b10^5`$ |
| LMSS,R=1/2, $`f`$=0.07 | 17.90 | 10.94 | 0.611 | 0.612 | $`p_b10^5,p_b10^5`$ |

## IV Complexity per iteration

In this section we show that the complexity per single iteration is almost the same for both methods. Hence, the gain in iterations represents the gain in decoding complexity.

For both schemes, the calculation of the pseudo posterior probabilities, the clipping and the convergence tests are identical, so they can be excluded from our discussion. Furthermore, the vertical passes in the two schemes are identical, hence the only remaining source for a possible difference in the complexity for the two schemes is the horizontal pass. For simplicity, we assume a regular matrix $`H`$ of dimensions $`m`$ rows by $`n`$ columns, which has $`k`$ nonzero elements per row and $`c`$ nonzero elements per column. (One can easily extend the discussion to include irregular matrices, but the conclusions are the same).

In the PUS, each horizontal pass consists of $`k`$ subtraction operations to find $`\delta q_{ij}`$’s and $`(k1)`$ multiplications to find $`\delta r_{ij}`$’s (using (6)) for each bit. The calculation of $`r_{ij}^1`$ & $`r_{ij}^0`$ from $`\delta r_{ij}`$ requires two additions and two multiplications (7). The total number of operations per iteration for the PUS is given by

$$m(k+2k)=3mk$$
(12)
$$m(k(k1)+2k)=mk(k+1)$$
(13)
In the SUS, horizontal passes are done separately for each bit, summing to $`nc`$ passes in total. Each horizontal pass consists of $`k1`$ subtractions to find $`\delta q_{ij^{}}`$ for all participants in the check, except for the current bit, and $`k1`$ multiplications are required to calculate $`\delta r_{ij}`$. The calculation of $`r_{ij}^1`$ & $`r_{ij}^0`$ from $`\delta r_{ij}`$ in this scheme requires two additions and two multiplications per bit. Hence the total complexity is given by

$$nc(k1+2)=nc(k+1)$$
(14)
$$nc(k1+2)=nc(k+1)$$
(15)
Recalling that $`mk=nc`$, we have:

* (13) equals (15), so the same number of multiplications is performed in both cases.
* (14) becomes $`m(k^2+k)`$, which for $`k>2`$ is larger than (12). However, for small $`k`$ ( $`k`$5 for KS, and $`k`$ is of the same order for LMSS), the increment in the total fraction of additions is small. Furthermore, one must remember that *the complexity is dominated by the number of multiplications*.

Note that the abovementioned comparison was made under a straightforward implementation of both algorithms. In advanced algorithms the following improvements can be adopted in order to reduce the complexity of the schemes.

1. Some savings can be made for the horizontal passes in the PUS, for instance, computing $`\delta q_{ij}`$ for the entire row, and dividing by each $`\delta q_{ij}`$ element, or recomputing $`\delta q_{ij}`$ only for updated bits in the SUS.
2. PUS can be implemented simultaneously over all checks (variables) using several processors. The implementation of SUS in parallel over a finite fraction of the checks (variables) is possible, but may require a special design.
3. SUS has some advantage in memory requirement, since only a column vector of the currently updated $`r_{ij}`$ is required, whereas for the PUS the whole $`r_{ij}`$ matrix must be retained simultaneously.

## V Qualitative theoretical explanation

The key difference between the two algorithms is the inter iteration information exchange, which is a property of the SUS only. Let us denote by $`r_{ij}^t,q_{ij}^t`$ the values computed in iteration $`t`$. In the PUS all $`r_{ij}^t`$ values are determined by the $`q_{ij}^{t1}`$ values (values of the previous iteration), and the $`q_{ij}^t`$’s are determined by these $`r_{ij}^t`$’s . In the SUS, after a bit is updated, the following bits that share a check with it are already exposed to the updated information. For instance, assume $`x_j`$ and $`x_k`$ share a check $`i`$, i.e. $`H_{ij}=H_{ik}=1`$, and assume $`j<k`$. In iteration $`t`$, $`r_{ij}^t`$ is updated using $`q_{ik}^{t1}`$; however, proceeding to the $`k`$th column, $`r_{ik}^t`$ is updated using $`q_{ij}^t`$, the most recent available information.

In other words, in the SUS, the first bits to be updated utilize the previous iteration data $`(q_{ij}^{t1})`$. A group of bits use mixed data from previous and current iterations and, finally, some of the bits are entirely updated by information from the current step $`(q_{ij}^t)`$.

The gain in the number of iterations for the SUS can be qualitatively well understood by the following argument: Since the decoding procedure terminates successfully (with some small $`p_b`$) and the number of correct bits increases monotonically, it is evident that on average, the current knowledge is superior to the knowledge of previous iteration.111Our simulations indicate that the number of correct bits versus time (iterations) is a *concave function* (except for the last few iterations). This result leads to a faster convergence for the SUS (many-small-steps converging strategy) in comparison to the PUS.

An important question is raised regarding which part of the SUS is accelerated in comparison to the PUS. The acceleration of the SUS may be a result of one of the following regimes: (a) a faster asymptotic convergence; (b) a faster arrangement in the initial stage of the decoding from random initial conditions; or (c) a uniform acceleration over all the stages of the decoding.

In order to answer this question we perform the following advanced simulations. We run the PUS and record the number of correct bits in each iteration. The correction gain of each iteration is defined as the increment in the fraction of correct bits. At each step of the PUS we prepare another replica of the system with the same initial conditions, the same $`q_{ij}`$ and $`r_{ij}`$, and run one iteration of the SUS. The correction gain of the SUS is then compared to that of the PUS. In Figure 3 we plot the rate between the sequential and parallel correction gains as a function of time (marked $`\times `$). This rate is nearly $`2`$, with relatively small fluctuations along the decoding process. In other words, on the average the SUS corrects twice the number of bits in comparison to the PUS, independent of the state of the decoder. (The simulation was performed on the KS construction with rate $`1/3`$, $`f=0.155`$, block size $`10000`$, $`20`$ different samples, and all convergence times were normalized to a 0-1 scale.)

The observation that the correction gain uniformly distributed over all the stages of the decoding raises the question of whether there is a superior updating order of the bits producing a correction gain greater than $`2`$. For one iteration of the KS construction, it may be better to update the bits from left to right than in the reverse order. For rate $`1/3`$, for instance, the right-most part of the matrix (about 25% of the columns) contain only one non-zero element per column, and this element is also the last element in its check. These bits are entirely updated by current iteration data, resulting in an increased correction gain. At the left-most end, on the other hand, there are 7 non-zero elements per column (for rate 1/3 construction), so that only a small fraction of them are the "last bit" for all the checks in which they participate. Most of the bits are updated by mixed information from the current and previous iteration, resulting in a smaller correction gain. In Figure 3 the rate between the correction gain for SUS and PUS for reversed (right to left) updating order is marked “$``$”. This rate is evidently less than for the left to right updating. Preliminary simulations indicate that by carefully selecting the updating order one can save 10%-30% of the iterations relative to a plain left to right sequential updating.

## VI Simulations

In this sections the technical details of our simulations are described. We generate the $`H`$ matrix at random, distributing the non-zero elements as evenly as possible without violating the constraint of the number of elements per row/column. No special attempt was made to select a "good performing" matrix. For the KS structure, we generate the $`x`$ vector as follows: The source bits were set to $`1`$ or $`0`$ with probability 0.5. The noise bits were set to $`0`$, and then exactly a fraction $`f`$ of the bits were selected randomly and flipped ($`f`$ is the flip probability). The check vector $`z`$ was computed by $`z=Hx`$, and the algorithm solved $`Hx^{}=z`$. We found $`p_b`$ by comparing $`x`$ & $`x^{}`$, for the source region only. The source length selected was $`n=10,000`$ (resulting in $`x`$ of length $`40,000`$, and $`z`$ of length $`30,000`$ for rate $`1/3`$).

For the LMSS structure, following shokrolahi spilman irreg 65-85  we always decoded the all-zero codeword, generating the noise vector $`n`$ in the same way as described above. The check vector $`z`$ was computed, $`z=Hn`$, and the algorithm solved $`Hn^{}=z`$. We found $`p_b`$ by comparing $`n`$ & $`n^{}`$ (in the LMSS version the "decoding" ends when the noise vector is found and the transmitted vector, $`t`$, is related to the received vector, $`r`$ by $`t=r+n`$(mod 2). Finding the source message from $`t`$ is not defined as part of the decoding problem). We used a noise vector of length $`20,000`$, corresponding to a check vector of size $`10,000`$ (rate 1/2).

In both cases the flip rate, $`f`$, was selected as being close enough to the critical rate for this block length such that the decoding is characterized by relatively long convergence times. However, the flip rate $`f`$ was chosen not too close to the threshold in order to avoid a large fraction of non-converging samples. After the check vector $`z`$ was constructed, it was decoded both in parallel and sequential schemes, and the number of iterations was monitored. We defined 3 halting criteria for the iterative process:

1. The outcome $`x^{}`$ fully solves $`Hx^{}=z`$.
2. The algorithm reached a stationary state, namely, $`x^{}`$ did not change over the last 10 iterations.
3. A predefined number of iteration was exceeded ("non-convergence"). This number was selected so as to be far larger than the average converging time (500 iterations in our case).

The vast majority of samples converged successfully. More precisely, less than 0.2% samples failed to converge or reach a non-solving stationary state.

## VII Conclusions

We demonstrated that the SUS outperforms the PUS in the convergence time aspect by a factor of about 2. Since the complexity per iteration of the two schemes is nearly the same, the gain in iterations is similar to the gain in the decoding complexity. The time gain is probably related to the inter-iteration information exchange, which is a property of the SUS. This explanation is also consistent with the observation that the gain is uniformly distributed over all the decoding stages. The question of whether the number of iterations can be reduced by a factor greater than $`2`$ by updating the bits in a special order is currently under investigation.

We acknowledge fruitful discussions with D. Ben-Eli and I. Sutskover.


text:
# QCD String Structure in Vector Confinement

## I Introduction

We present here an interesting simple observation on electric potential confinement. Our conclusions concern the nature of potential models and their relationship to QCD string confinement. In the case of the time component vector (TCV) potential we note a remarkable transition in the Regge spectroscopy from the leading trajectories with slopes $`1/4a`$ to deep daughter slopes of $`1/\pi a`$. In fact, the exact string spectroscopy is found to be embedded in this regime. The latter slope is characteristic of the Nambu-Goto QCD string with one fixed end ref:one . We first discuss some general properties of the spinless Salpeter (SS) or square root equation. For the linear TCV potential we then compute the Regge structure in the orbitally and radially dominant regimes.

## II Spinless Salpeter Equation

The spinless Salpeter equation has long been attractive as a minimal generalization of the non-relativistic Schrödinger equations with relativistic kinematics. Its main difficulty is that it involves the square root of a differential operator, sometimes called a pseudo-differential operator, which is difficult to treat analytically. Even so, it is easily implemented numerically. The introduction of a variational method ref:two  allowed straightforward numerical solutions involving square roots and even transcendental functions of operators. In this section we discuss, though do not prove, the necessity of using the spinless Salpeter equation in place of a Klein-Gordon type equation when dealing with a pure vector potential.

Starting from the simple classical Lagrangian,

$$L=(m+S(r))\sqrt{1v^2}V(r),$$
(1)
describing a particle of mass $`m`$ moving with velocity $`v`$ in a Lorentz scalar potential $`S(r)`$ and a Lorentz time-component vector potential $`V(r)`$, we find that the Hamiltonian,

$$H=\pm \sqrt{p^2+(m+S(r))^2}+V(r),$$
(2)
has a sign ambiguity. This ambiguity arises from the elimination of the velocity in favor of the momentum

$$𝐩=\frac{L}{𝐯}=(m+S(r))\gamma 𝐯,$$
(3)
where $`\gamma =(1v^2)^{1/2}`$. The positive sign in Eq. (2) corresponds to the usual positive energy states and yields the correct non-relativistic limit

$$H\stackrel{\left|p\right|m}{}m+\frac{p^2}{2m}+S(r)+V(r).$$
(4)
The squared form of Eq. (2),

$$(HV(r))^2=p^2+(m+S(r))^2,$$
(5)
known as the generalized Klein-Gordon (GKG) equation, behaves quite differently in the cases of a pure vector and a pure scalar potential. This difference is most easily seen from consideration of the s-wave turning points, where $`p^2=0`$. Eq. (2) determines the energy as a function of the potentials at the turning points,

$$E=\pm (m+S)+V.$$
(6)
We show this turning point condition separately for pure scalar and pure vector linear confining potentials in Figs. 1 and 2. In Fig. 1 we see that for the pure scalar case there is no overlap of the positive and negative branches of Eq. (6). In Fig. 2, on the other hand, we see that for any $`E>m`$ in pure TCV linear confinement there are states on both positive and negative branches. A completely analogous observation ref:three  has been made earlier in the context of the Dirac equation.

Since the GKG equation (5) is quadratic and contains both branches of Eq. (2), it can have no normalizable solutions in pure TCV confinement because there is tunneling to the negative energy states, which becomes catastrophic for small quark mass. This difficulty is the well-known “Klein Paradox” ref:four . On the other hand, the GKG equation with pure scalar confinement does have well defined normalizable solutions because there is no possibility of tunneling to the negative branch.

The spirit of the Salpeter equation ref:five  is to conserve a definite particle number. To this end, energy projection operators are employed to remove the negative energy spectrum and the fermion TCV Salpeter equation has normalizable solutions for confined quarks ref:three . The spinless Salpeter equation for TCV confinement achieves the same end more simply by explicitly including only the positive branch of Eq. (2). For p-wave and higher angular momentum states, we consider the Hamiltonian in spherical coordinates,

$$H=\sqrt{p_r^2+\frac{J^2}{r^2}+m^2}+ar,$$
(7)
which we obtain from dividing the momentum into radial and angular parts. The wave equation to be solved is then the eigenvalue problem

$$H|\psi =E|\psi $$
(8)
for the pseudo-differential operator given by Eq. (7) with the usual operator replacements.

## III Semi-Classical Analytic Quantization

We have discussed the necessity of using the SS equation instead of the GKG equation with vector confinement. In semi-classical language, the necessity of using the SS equation becomes the necessity of using the turning points given by Eq. (6) with the positive sign only. In this section we calculate the semi-classical spectrum of the SS equation with TCV confinement. We first rewrite the SS equation (7) in terms of dimensionless variables

$`x`$ $``$ $`{\displaystyle \frac{ar}{E}},`$
$`\beta `$ $``$ $`{\displaystyle \frac{aJ}{E^2}},`$ (9)
$`P_r`$ $``$ $`{\displaystyle \frac{p_r}{E}}.`$
In the case of a massless quark, we have

$`P_r`$ $`=`$ $`{\displaystyle \frac{\sqrt{Q}}{x}},`$ (10)
$`Q`$ $``$ $`(xx_1)(xx_2)(x_{1+}x)(x_{2+}x),`$
where

$`x_{1\pm }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1\pm \sqrt{14\beta }\right),`$
$`x_{2\pm }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(1\pm \sqrt{1+4\beta }\right).`$ (11)
The dimensionless physical turning points resulting from the SS equation are $`x_{1\pm }`$. The dimensionless constants $`x_{2\pm }`$ are turning points of the GKG equation, but not of the SS equation.

In the regime we are investigating the angular momentum $`J`$ is limited but the energy $`E`$ is large. It follows that $`\beta 1`$. To see this, we note that even for the leading (circular orbit) Regge trajectory $`J/E^2=1/4a`$ (see the appendix) and hence $`\beta _{\mathrm{leading}}=1/4`$. For fixed $`J`$ and large $`E`$, $`\beta `$ becomes small. In this limit the leading behavior of the roots (III) is

|  | $`x_1=\beta +\beta ^2+\mathrm{},`$ | $`x_{1+}=1\beta \beta ^2\mathrm{},`$ |
| --- | --- | --- |
|  | $`x_2=\beta +\beta ^2+\mathrm{},`$ | $`x_{2+}=1+\beta \beta ^2+\mathrm{}.`$ |

(12)
The SS equation may be quantized semi-classically by the usual method ref:six

$$_r_{}^{r_+}𝑑rp_r=\pi \left(n+\frac{1}{2}\right).$$
(13)
In dimensionless variables (III), the quantization condition (13) becomes

$$_x_{}^{x_+}𝑑xP_r(x)=\frac{a\pi }{E^2}\left(n+\frac{1}{2}\right).$$
(14)
The TCV quantization is achieved by integrating Eq. (10). The method is quite accurate for all states and becomes exact for states with many nodes in the radial wavefunction. This is the regime we are particularly interested in.

Although the integral of Eqs. (10), (14) can be expressed in terms of elliptic integrals, it is more efficient to approximate the integrand first for small $`\beta `$. The approximation

$$\frac{\sqrt{Q}}{x}\frac{\sqrt{(xx_1)(xx_2)}}{x}\left[1x\frac{\beta ^2}{2(1x)}\beta ^2\right]$$
(15)
reproduces the exact integrand to better than $`0.5\%`$ throughout most of the region of integration, even when $`\beta =0.1`$. The accuracy of this approximation is shown graphically in Fig. 3, where we plot both the exact integrand (10) and the approximation (15) for $`\beta =0.1`$.

With the approximation (15), the quantization integral (14) is relatively easily evaluated. To order $`\beta `$ we find

$$_x_{}^{x_+}\frac{\sqrt{Q}}{x}𝑑x\frac{1}{2}\frac{\pi \beta }{2}+𝒪(\beta ^2),$$
(16)
which immediately leads to the result

$$\frac{E^2}{\pi a}=J+2n+1.$$
(17)
Upon making the Langer ref:seven  correction to take into account the centrifugal singularity,

$$JJ+\frac{1}{2},$$
(18)
where $`J`$ is now the angular momentum quantum number, we find the final spectroscopic relation

$$\frac{E^2}{\pi a}=J+2n+\frac{3}{2}.$$
(19)
The relation (19) is identical to an analytic approximation for highly radially excited QCD string, as well as being an excellent approximation to its exact numerical solution ref:eight .

## IV Comparison to Exact Numerical Results

We can nail down our central result that the QCD string spectrum (19) is replicated in the radially dominant regime of the TCV potential by solving the TCV SS equation exactly numerically. The variational (Galerkin) method is well suited for solving eigenvalue equations with mixed coordinate and momentum operators. Briefly, the method begins with a complete set of orthogonal states that can be Fourier transformed. The wavefunction is approximated as a superposition of the lowest $`N`$ of those states. The wave equation transforms into an $`N\times N`$ matrix equation which is then diagonalized. The accuracy of the resulting eigenvalues and wavefunctions is measured by the dependence on $`N`$ and the dependence of the scale parameter of the basis set. Details can be found in ref:nine . Some further considerations of the accuracy of the variational method are presented in ref:ten .

The TCV SS equation (7) with $`m=0`$ can thus be solved for a variety of angular momenta, $`J=0,1,2,\mathrm{}`$ and radial states, $`n=0,1,2,\mathrm{}`$. The result is depicted in Fig. 4. Plotted on this figure are both the QCD string slope $`1/\pi a`$ (solid) and the usual leading trajectory TCV slope $`1/4a`$ (dashed). In both cases we start the lines at $`J=0`$. Although the leading trajectory ($`n=0`$) agrees with the dashed prediction, the situation changes as we examine the higher radial excitations. At the larger $`n`$ values (deep daughters) we observe that the solid line corresponding to the QCD string slope becomes accurate.

In Table 1 we provide the exact and the string solutions for the s-wave radial excitations. We see now that all aspects of the QCD string result (19) work quite well. The difference between excitations is accurately two units of $`E^2/\pi a`$ and the absolute values of the state energies correspond well to the QCD string to three significant figures.

## V Summary and Discussion

On occasion, some simple results are quite unexpected. We did not expect to find the $`m=0`$ QCD string spectroscopy in a potential model. To summarize our finding, we have examined both numerically and analytically the state spectroscopy for a massless quark moving in the time-component vector linear confinement potential. For states with angular momentum much less than the number of nodes we find the spectrum is exactly that of the QCD string. The Regge slope, radial excitation energy, and the absolute values of the energy (19) are exactly what one expects from the QCD string.

Simple results, even surprising ones, usually have simple explanations. As we have recently pointed out ref:eight , the QCD string equations reduce to the spinless Salpeter equation with a linear time-component vector potential for the s-wave states. Thus there is a natural physical connection between the two systems. We can understand this connection simply. A QCD string corresponds to a constant chromoelectric field in the quark rest frame. In the limit that the quark is moving radially, the QCD string has no angular momentum. The total energy of the system becomes

$$E_{\mathrm{string}}=\sqrt{𝐩^2+m^2}+ar\frac{\mathrm{arcsin}v_{}}{v_{}}\stackrel{v_{}0}{}\sqrt{p_{r}^{}{}_{}{}^{2}+m^2}+ar.$$
(20)
In a potential model the field never carries angular momentum. The potential energy of the quark is linear in the distance from the origin, which is the same as the energy of the string as long as the string is moving radially so that there are no relativistic corrections due to its transverse motion. The total energy of the quark in a TCV potential,

$$E_{\mathrm{TCV}}=\sqrt{𝐩^2+m^2}+ar,$$
(21)
in the limit of vanishing $`v_{}`$ becomes

$$\underset{v_{}0}{lim}E_{\mathrm{string}}=\underset{v_{}0}{lim}E_{\mathrm{TCV}}=\sqrt{p_{r}^{}{}_{}{}^{2}+m^2}+ar.$$
(22)
It is thus natural to expect agreement between string and TCV confinement in the radial dominant regime.

The real puzzle is why the Regge slope of the QCD string should be the same for both circular and radial motions. On one hand it is well known from Nambu-Goto days that the Regge slope of a rotating QCD string with one end fixed is $`1/\pi a`$. We have seen that the TCV potential Regge slope varies from $`1/4a`$ to $`1/\pi a`$ as one goes from orbital to radial motion \[see Fig. 4\]. There is thus no obvious reason why this should not happen with the QCD string. Since the QCD string/TCV potential Regge slopes coincide at $`1/\pi a`$ for nearly radial motions this explains the remarkable uniform Regge structure seen previously from numerical solutions of the QCD string ref:eight .

## Acknowledgments

This work was supported in part by the US Department of Energy under Contract No. DE-FG02-95ER40896.

## Appendix A Nearly Circular Orbits

For completeness we semi-classically quantize the TCV equation for nearly circular orbits ($`Jn`$). Using the notation of Eqs. (IIIIII), we consider the orbitally dominant regime. In this case the turning points $`x_{1+}`$ and $`x_1`$ are nearly equal. By Eq. (III), this occurs exactly at

$`x_c`$ $`=`$ $`{\displaystyle \frac{1}{2}},`$
$`\beta _c`$ $`=`$ $`{\displaystyle \frac{1}{4}}.`$ (23)
For nearly circular orbits we can expand in $`\beta `$ and $`x`$, giving

$`x_{1\pm }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\pm \sqrt{{\displaystyle \frac{1}{4}}\beta }`$
$`Q`$ $``$ $`{\displaystyle \frac{1}{16}}\beta ^2{\displaystyle \frac{1}{2}}\left(x{\displaystyle \frac{1}{2}}\right)^2`$
$``$ $`{\displaystyle \frac{1}{2}}(xx_1)(x_{1+}x)`$
Using the quantization condition (14), we obtain

$$\frac{1}{2\sqrt{2}}\left(12\sqrt{\beta }\right)=\frac{a}{E^2}\left(n+\frac{1}{2}\right),$$
(25)
with $`\beta =aJ/E^2`$. Solving for $`\sqrt{\beta }`$, squaring and dropping the small squared $`n/E^2`$ term, we obtain

$$\frac{E^2}{4a}=J+\sqrt{2}n+\frac{1}{\sqrt{2}}.$$
(26)
Finally, we make the Langer ref:seven  correction $`JJ+\frac{1}{2}`$ and find ref:eight

$$\frac{E^2}{4a}=J+\sqrt{2}n+\frac{1}{2}+\frac{1}{\sqrt{2}}.$$
(27)