text:
# Commensurate-incommensurate transitions in quantum films: Submonolayer molecular hydrogen on graphite

## I Introduction

The commensurate-incommensurate transition has been expensively studied in classical monolayer films both theoretically and experimentally. Monolayers of hydrogen or helium are quantum mechanical systems and, in principle, one might suspect a different behavior due to the effect of strong quantum fluctuations.

The monolayer phase diagram of molecular para-hydrogen adsorbed on graphite as inferred from experimental studiesFreimuth87  is shown in Fig. 1. This phase diagram was originally drawn based on the anomalies found in the specific heatMotteler ; Chaves ; Freimuth85 ; for the characterization of the various phases, low energy electron diffraction (LEED)Cui86 ; Cui  as well as neutron diffractionNielsen77 ; Freimuth87 ; Wiechert91a  studies have been carried out. At 1/3 coverage the molecules condense on the surface of graphite in a $`\sqrt{3}\times \sqrt{3}`$ commensurate solid. In the low coverage region $`(\rho <0.6)`$ (in units of the $`\sqrt{3}\times \sqrt{3}`$ commensurate density), it forms a commensurate solid-gas coexistence phase at low temperature and a 2D gas phase at higher temperature ($`T>10K`$). For coverages $`0.6_{}^<\rho _{}^<0.9`$, as a function of temperature, there is a transition from a commensurate solid cluster phase to a commensurate solid phase with vacancies at higher temperature and to a 2D gas at even higher temperature. At density somewhat higher than unity the so-called $`\alpha `$-phase is formed which is believed to be a striped domain-wall solid phase and at higher temperature it transforms to the so-called $`\beta `$-phase and to a fluid phase at even higher temperature. At high densities a transition take place to a triangular incommensurate solid phase.

The phase diagram of molecules formed from the isotopes of the hydrogen such $`D_2`$ and $`HD`$ molecules physiorbed on graphite has also been studiedFreimuth90 ; Lauter90  and it is similar to that of $`H_2`$ on graphite but more complex.

We have recently studiedNM  the phase diagram of molecular hydrogen on graphite at and below 1/3 coverage using Path Integral Monte Carlo (PIMC) simulation. Our computed phase diagram was in general agreement with that inferred from the experimental studies for this coverage range. In this paper we extend the PIMC calculation to study the phase diagram above the commensurate density which is less well understood both theoretically and experimentally. At coverages above $`\rho 1.05`$ the monolayer undergoes a commensurate-incommensurate(C-IC) transition. While the existence of the so-called $`\alpha `$ and $`\beta `$ phases were known from the specific heat anomalies, their characterization came from LEEDCui86 ; Cui  and neutron scatteringFreimuth87  studies. Using LEED experiments Cui and FainCui86 ; Cui  observed a uniaxial IC solid with striped superheavy domain walls and a rotated triangular IC solid at higher coverages $`\rho 1.23`$. Freimuth, Wiechert, and LauterFreimuth87  (FWL) presented a neutron diffraction study of the C-IC transition and their results are in agreement with the LEED results. They examined the commensurate-incommensurate solid transition, especially the striped domain-wall solid phase. In the $`\alpha `$-phase the diffraction intensity has a main peak characteristic of a compressed lattice and a satellite peak on the lower side of the commensurate peak wave-vector position $`k_c=1.703\AA ^1`$ that arise from the spatial modulation due to ordered striped-domain walls. As coverage increases, the separation of the two peaks increases and the satellite intensity decreases. At higher temperatures the peak height drops and the satellite peak vanishes which implies that the commensurate domains vanish and the system becomes an isotropic fluid phase. As coverage increases, the first molecular hydrogen layer forms a rotated incommensurate solid (RIC) phase. Neutron diffraction studies show a sharp and intense peak at wave-vector $`k=1.97\AA ^1`$ at the density $`\rho =1.34`$ which reveals an IC equilateral triangular structure. This phase is continuously compressed as the coverage increases up to the highest density allowed before layer promotion.

In parallel to these experimental investigations several important theoretical studies were undertaken to understand the commensurate-incommensurate transitionnovaco-mctague ; vilain ; bak ; Copper ; Timothy  in physiorbed systems. The phases and the phase transitions which were theoretically predicted and studied can be understood as classical phenomena due to the strain induced on the adsorbate molecules of the monolayer film due to the substrate periodic structure which arise from the competition between adsorbate-adsorbate and adsorbate-substrate interactions. Molecular hydrogen and atomic helium physisorbed on graphite are quantum films characterized by strong zero point motion. Therefore, one could question the degree to which a classical picture might be valid and might expect new phenomena and phases to occur. Motivated by these thoughts we extended our earlier investigationNM  to study this system above 1/3 coverage up to layer completion starting from the known hydrogen molecule-molecule and hygrogen-graphite interactionsSilvera ; Crowell ; Steele ; Carlos  and using the PIMCCeperley95  which is a Quantum Monte Carlo technique adopted for the study of strongly interacting quantum films by Pierce and ManousakisPierce .

We have simulated the first layer of molecular hydrogen on graphite with a variety of simulation cells that are appropriate for examining different phase regions of the phase diagram. We have applied periodic boundary conditions along both directions parallel to the surface of the substrate. We have computed expectation values for the total energy, the static structure factor, the probability distribution, and the specific heat by means of the path-integral Monte Carlo simulation method using the multi-level metropolis method. In the multi-level bisection method we have used level 3 which optimizes the acceptance ratio. To thermalize the system, we typically carried out of the order of 15,000 MC steps and we carried out of the order of 2,000-20,000 MC steps in order to compute observables.

## II Uniaxially Compressed Solid With Ordered Domain Walls

In this section we discuss the region above the commensurate coverage where experimentally the ordered domain-wall solid was found. We carried out a simulation at two densities in this region, the same densities discussed by FWL, where they inferred the existence of relaxed superheavy domain walls.

We used simulation cells that can accommodate the structures that have been proposed by FWL. a) First, we consider $`\rho =0.0694\AA ^2`$. We have chosen a simulation cell with dimensions $`x`$= 5$`\sqrt{3}`$ $`a_{gr}`$ and $`y`$= 22 $`a_{gr}`$, where $`a_{gr}=2.459\AA `$ is the carbon-carbon distance on the graphite surface and 80 hydrogen molecules to produce the above density. In Fig. 2 we give the contour plot of the calculated probability distribution where we see that two commensurate-solid domains are separated by the incommensurate solid domains. The solid structure is uniaxially compressed along our $`y`$ direction such that a new row of molecules is introduced for every 8 rows. Notice that the period in the $`x`$-direction is $`\sqrt{3}a_{gr}`$ and there is modulation along the $`y`$ direction with period $`11a_{gr}`$, so that the wavelength $`\lambda _s`$ of this striped domain-wall modulation is $`\lambda _s=27.049`$. b) Second, we have performed the simulation at the density $`\rho =0.0716\AA ^2`$. We have chosen a simulation cell with dimensions $`x=5\sqrt{3}a_{gr}`$ and $`y=16a_{gr}`$ and 60 hydrogen molecules in order to achieve the above mentioned density. The calculated contour plot of the probability distribution is shown in Fig. 3. Notice that in this case also there is an ordered striped domain-wall solid phase along our $`y`$ direction. Namely, the amplitude of the molecular density wave is modulated along the $`y`$ axis with a period of about half our cell size along the $`y`$ axis. Notice that there are two commensurate solid domains separated by denser regions.

Our computed static structure factor $`S(k)`$ for $`\rho =0.0694\AA ^2`$ is shown in Fig. 4. The main Bragg peaks of this structure occur at $`𝐤_1`$ = (1.475 $`\AA ^1`$,0.929 $`\AA ^1`$) and $`𝐤_2`$ = (0,1.858 $`\AA ^1`$). The satellite peak occurs at $`𝐤_{sat}^c`$= (0,1.626 $`\AA ^1`$). The experimental values of the magnitude of the wave-vectors at the peaks are $`k_{main}^{exp}`$= 1.743 $`\AA ^1`$ and $`k_{sat}^{exp}`$= 1.632 $`\AA ^1`$, which compare well with our computed values of $`k_1^c`$= 1.743 $`\AA ^1`$ and $`k_{sat}^c`$= 1.626 $`\AA ^1`$. We believe that FWL could not observe the peak at $`𝐤_2`$ because of the strong interference with the (002) graphite peak.

The interpretation of these results is as follows: Analyzing the contour plot of Fig. 2 we find that a) the solid appears uniaxially compressed along the $`y`$ direction, by adding another row for every 10 rows of molecules and spreading them evenly, while along the $`x`$ direction the average spacing between the molecules remains that of the commensurate solid. b) superimposed there is a density modulation along the $`y`$ direction which has wavelength several times greater than the average nearest neighbor distance but of small amplitude. The actual size of the wavelength $`\lambda _s`$ of this striped domain-wall modulation is half of the length of our simulation cell in the $`y`$-direction, i.e., $`\lambda _s=27.05\AA `$, as discussed in the previous paragraph. The two main Bragg peaks correspond to the unit vectors which span the reciprocal space of this uniaxially compressed triangular solid structure. The satellite peaks are due to lateral modulation in our $`y`$ direction and are located at $`𝐤_{sat}`$ = $`𝐤_{1,2}𝐤_s`$. $`𝐤_{1,2}`$ correspond to the wave vector of each of the two main Bragg peaks and $`𝐤_s=(0,2\pi /\lambda _s)`$, where $`\lambda _s`$ is the wave length of a unit cell for the striped domain-wall solid structure in the modulated direction. Using this $`\lambda _s`$ value, $`𝐤_2=(0,1.858\AA ^1)`$ and $`𝐤_s=(0,0.232\AA ^1)`$ (obtained by using the values of $`\lambda _s=27.049\AA `$ found by inspection of Fig. 2) we can find that the satellite peak position is at $`𝐤_{sat}=(0,1.626\AA ^1)`$, which agrees well with our peak value.

The structure factor for the case of $`\rho =0.0716\AA ^2`$ (Fig. 3) is shown in Fig. 5. The main peaks are $`𝐤_1`$ = (1.475 $`\AA ^1`$,0.958 $`\AA ^1`$) and $`𝐤_2`$ = (0,1.916 $`\AA ^1`$) which correspond to the uniaxially compressed solid produced by adding another row of molecules for every 6 rows of molecules and spreading them evenly. This, in addition, forms a superscructure with wavelength $`\lambda _s=19.672\AA `$ which gives a satellite peak at $`𝐤_{sat}^c`$= (0,1.597 $`\AA ^1`$). Notice that the scattering intensity, relative to the previously discussed case of $`\rho =0.0694\AA ^2`$, has the peaks shifted at higher values of $`k`$ as expected and the intensity of the satellite peaks is decreased as in the experiments.

## III Domain Wall Melting

In Fig. 6 we show the contour plot of the the probability distribution where a striped domain-wall fluid phase (at $`T=11.11K`$) and a fluid phase (at $`T=16.67`$) are evident. Notice that the stripe domain-wall fluid phase (left part of Fig. 6) is characterized by mobile commensurate and incommensurate domains. To understand this contour plot we need to be reminded of the periodic boundary conditions used in our simulation and the fact that these domains cannot intersect each other. This implies that one domain will oscillate back and forth between the neighboring domains. We can determine the melting temperature of the domain-wall solid phase from the temperature dependence of the specific heat and from the temperature dependence of the peak height of the static structure factor $`S(k)`$.

The calculated specific heat as a function of temperature is shown in Fig. 7 and is characterized by two peaks, the first corresponds to the melting of the domain-wall solid while the second peak indicates the melting of the stripes and the solid into a fluid. Our computed values for the melting temperature are generally lower that the experimental values. For example, for $`\rho =0.0716\AA ^2`$ we find that at $`T=9.09K`$ the striped domain-wall solid phase undergoes a transition to the domain-wall fluid phase and at $`T=12.5K`$ the $`\beta `$-phase becomes an isotropic fluid phase. Factors for obtaining lower values for the critical temperature than the experimental values could be the finite-size effects and the interaction strength used to describe the hydrogen-graphite interaction. We also computed the static structure factor shown in Fig. 8 and we studied the temperature dependence of the height of its first main peak which is shown in Fig. 9. The peak height significantly decreases near the melting temperature of the domain-wall solid and near the domain-wall evaporation temperature.

## IV Rotated Incommensurate Solid

The first layer of molecular hydrogen adsorbed on graphite forms an incommensurate solid phase before the first layer is complete. We have carried out a simulation at the density $`\rho =0.0849\AA ^2`$ using a simulation cell with dimensions $`x`$= 9 $`\sqrt{3}`$ $`a_{gr}`$ and $`y`$= 9 $`a_{gr}`$. The calculated probability distribution(Fig. 10) also clearly shows the presence of an equilateral triangular solid structure which is not registered with the underlying graphite lattice and it is rotated relative to the graphite lattice. This rotation was first predicted and discussed by Novaco and McTaguenovaco-mctague . The angle of the incommensurate lattice relative to the graphite lattice is approximately $`5^{}`$ in good agreement with the experimental valueCui  for the above density. The calculated static structure factor (shown in Fig. 11) has one sharp peak at $`k`$= 1.99 $`\AA ^1`$ corresponding to an unregistered equilateral triangular lattice in agreement with the value of $`k`$= 1.97 $`\AA ^1`$ reported from the neutron diffraction experimental studyFreimuth87 .

## V Energy Calculation

In Fig. 12 we compare the energy per hydrogen molecule as a function of the surface coverage for the cases where the full potential is used (solid line) and the case where only the laterally averaged potential is used (dashed line). There are several important observations which will make based on the results presented on the this figure. First the curve obtained using the laterally averaged potential has a minimum at the density $`\rho _0=0.0705\AA ^2`$. At this equilibrium density the system forms a triangular lattice which fills the entire system and below this density the system undergoes a phase separation and forms a solid-vapor coexistence phase. Therefore below $`\rho _0`$ the energy as a function of density is higher for a finite-size system or for a system which is forced to be uniform at this lower density. Above this density the system is a compressed triangular solid until the promotion densityNM . The solid line corresponds to the energy per molecule when we used the hydrogen-molecule-graphite interaction potential with the corrugations. The minimum in this case occurs at $`\rho _c=0.0636\AA ^2`$ which corresponds to the $`\sqrt{3}\times \sqrt{3}`$ solid. Below that density the system is unstable to formation of solid clusters and this part of the phase diagram was investigated in Ref.NM . Above this density up to a density one can notice that the presence of the substrate corrugations moves the minimum from $`\rho _0`$ at much higher energy but there is a feature at the same density produced by the $`H_2`$-$`H_2`$ interaction. The value $`\rho _0`$ of the equilibrium density for a smooth substrate plays a significant role even when the full potential with the substrate corrugations is used because it is the density which minimizes the energy due to $`H_2`$ molecule-molecule interaction. Clearly due to this term which has a preference for a higher density than the commensurate density, the solid line is not convex which a sign of instability of a uniform phase. These results are obtained on lattices which frustrate the proper periodicity of the domain walls. If the calculations are done on lattice which can accommodate the periodicity of these superlattice structure the energy is lowered. As an example we also give (solid circle) the energy obtained for the stripe state of Fig. 3. Notice that the energy falls below that obtained with choice of lattices unfavorable for this ordered domain-wall state.

## VI Conclusions

Our studies of $`H_2`$ monolayer on graphite above 1/3 coverage clarifies the nature of the phase diagram. The boundaries of the phase diagram were well-known from the specific anomalies. We find that the phase near 1/3 coverage is the known $`\sqrt{3}\times \sqrt{3}`$ commensurate solid. At higher densities the monolayer undergoes a transition to a uniaxially compressed solid with density wave stripe-like modulations along the same direction. This phase was first identified by LEEDCui86  studies; while the subsequent neutron diffraction studiesFreimuth87  confirm the existence of ordered domains walls, the peak which is the signature of the uniaxial compression is not seen because of the (002) reflection from the underlying graphite lattice. In addition in both studies the domain walls are assumed to be superheavy. As it is clear from both our contour plots and our structure factor the walls are not significantly denser than the rest of the system; we find that the additional molecules added to the monolayer to increase the coverage above the commensurate coverage are used primarily to uniaxially compressed the solid and they are only partially used to construct the walls. Our calculated specific heat for this coverage range as a function of temperature has two anomalies at two temperature values. From the calculated contour plots and the structure factor, the lower temperature anomaly is identified as the transition from the ordered domain-wall solid to a solid where the domain walls form a fluid and the higher temperature anomaly is caused by the total melting of the solid to a uniform fluid.

At even higher density we find that the uniaxially compressed solid with the ordered domain walls undergoes a transition to a triangular solid which is incommensurate with respect to the substrate lattice and in addition the triangular solid of the adlayer is rotated by an angle of the order of $`5^{}`$ with respect to the underlying substrate lattice. This is in complete agreement with both the theoretical predictionnovaco-mctague  and experimental findingsCui86 ; Freimuth87 .

## VII Acknowledgements

This work was supported by a National Aeronautics and Space Administration grant number NAG8-1773.


format: text

text:
# A strong astrophysical constraint on the violation of special relativity by quantum gravity11footnote 1Previous title: Lorentz violation and Crab synchrotron emission: a new constraint far beyond the Planck scale

## Methods

To determine the effects of Lorentz violation of (3) we follow the heuristic derivation of (3) given in Ref.  assuming the framework of effective field theory. Without assuming Lorentz invariance, the purely kinematical result is

$$\omega _c=\frac{3}{4}\frac{1}{R(E)\delta (E)}\frac{1}{c(\omega _c)v(E)}$$
(8)
where $`R(E)`$ is the radius of curvature of the orbit, $`\delta (E)`$ is the opening angle for the forward-directed radiation pattern, and $`c(\omega _c)`$ and $`v(E)`$ are the group velocities of the radiation and electron respectively. The solution of Eq. (8) for $`\omega _c(E)`$ determines the cutoff synchrotron frequency. In Eq. (8) we have used the fact that the electron and photon speeds are very close to the low energy speed of light $`c`$, which is set equal to unity. The numerical constant is chosen to yield the correct relativistic result (3) with $`\delta (E)=\gamma ^1(E)`$ (where $`\gamma (E)=\gamma (v(E))`$).

The radius $`R(E)`$ for a given energy is determined by the equation of motion of the electron in a magnetic field. All Lorentz violating terms in the equation of motion are suppressed by $`M^1`$, and the leading high energy corrections come from modifications to the minimal coupling terms. To estimate the magnitude of the change we use the dispersion relation as the Hamiltonian. As usual the minimal coupling is incorporated replacing the momentum by $`𝐩e𝐀`$, where $`𝐀`$ is a vector potential for the magnetic field. This yields the equation of motion $`𝐚=[1+3\eta E/2M](e/E)𝐯\times 𝐁`$, where we have kept only the lowest order term in $`\eta `$ and assumed relativistic energy $`Em`$. Since $`EM`$, the presence of the Lorentz violation makes very little difference to the orbital equation, hence we conclude that to a very good approximation the radius is related to the magnetic field and the energy of the electron by the standard formula $`R(E)=E/eB`$ (where again the speed of the electron has been set equal to unity).

The angle $`\delta (E)`$ scales in the Lorentz invariant case as $`\gamma ^1(E)`$. Since $`R(E)`$ and hence the charge current is nearly unaffected by the Lorentz violation, any significant deviation in $`\delta (E)`$ can only come from the modified response of the electromagnetic field to a given current. The amplitude of the modification is suppressed by at least one power of $`M`$ however, so dimensional analysis implies that it is suppressed by something of order $`\xi \omega /M`$. For the Crab synchrotron radiation used in our constraint $`\omega `$ is 100 MeV, which yields a suppression factor of $`10^{20}\xi `$. Therefore the amplitude of the electromagnetic field modification is negligible in comparison to the zeroth order synchrotron radiation (which is still present), and hence $`\delta (E)`$ scales with $`\gamma ^1(E)`$ in the usual way.

Since the emitted photons have relatively low energy compared to the electrons, it turns out that $`\xi `$ can be neglected in the relevant region of parameter space. Thus, since $`v(E)`$ is very close to $`c`$, the reciprocal of the difference of group velocities in the last term of Eq. (8) is well approximated by $`2\gamma ^2(E)`$. This yields (4).

The maximum synchrotron frequency $`\omega _c^{\mathrm{max}}`$ is obtained by maximizing $`\omega _c`$ (4) with respect to the electron energy, which amounts to maximizing $`\gamma ^3(E)/E`$. The difference of group velocities is given by

$$c(\omega )v(E)=\xi \frac{\omega }{M}+\frac{m^2}{2E^2}\eta \frac{E}{M},$$
(9)
where we have used the dispersion relations (1,2) and dropped higher order terms. Dropping the $`\xi `$ term this yields

$$\gamma (E)\left(\frac{m^2}{E^2}2\eta \frac{E}{M}\right)^{1/2},$$
(10)
Using this expression we find

$$\omega _c^{\mathrm{max}}=0.34\frac{eB}{m}(\eta m/M)^{2/3},$$
(11)
which is attained at the energy $`E_{\mathrm{max}}=(2m^2M/5\eta )^{1/3}=10(\eta )^{1/3}`$ TeV. Since $`\eta `$ is negative this is higher than the Lorentz invariant value that produces the same frequency, but only by a factor of order unity that works out to be $`(9/5)^{3/4}1.55`$. The energy of the electrons that produce the synchrotron radiation of frequency 100 MeV is 1500 TeV in the Lorentz invariant case. If this were the maximum frequency in the Lorentz violating case, it would be produced by $`2300`$ TeV electrons.

We now analyze the question of neglecting $`\xi `$ in evaluating the difference of group velocities (9). At the energy $`E_{\mathrm{max}}`$, the ratio of the $`\xi `$-dependent term to the other terms is given by

$$\frac{|\xi |\omega }{(m^2M/2E^2)\eta E}=3\times 10^{11}|\xi |(\eta )^{4/3}.$$
(12)
Hence neglecting $`\xi `$ is justified provided

$$|\xi |10^{11}(\eta )^{4/3}.$$
(13)
## Acknowledgements

We wish to thank F.A. Aharonian, G.E. Allen, G. Amelino-Camelia, and F. Stecker for helpful discussions. This work was supported in part by the NSF.


text:
# On geometry effects in Rayleigh-Bénard convection

## Abstract

Various recent experiments hint at a geometry dependence of scaling relations in Rayleigh-Bénard convection. Aspect ratio and shape dependences have been found. In this paper a mechanism is offered which can account for such dependences. It is based on Prandtl’s theory for laminar boundary layers and on the conservation of volume flux of the large scale wind. The mechanism implies the possibility of different thicknesses of the kinetic boundary layers at the sidewalls and the top/bottom plates, just as experimentally found by Qiu and Xia (Phys. Rev. E58, 486 (1998)), and also different $`Ra`$-scaling of the wind measured over the plates and at the sidewalls. In the second part of the paper a scaling argument for the velocity and temperature fluctuations in the bulk is developed.

Turbulent Rayleigh-Bénard convection is one of the classical problems in fluid dynamics . The great interest into this problem presumably also originates from the relevance of thermal turbulence in meteorology, geophysics, oceanography, and astrophysics. However, in contrast to these natural problems, in laboratory experiments the thermal convection is confined to a container. Recent experiments have revealed that geometrical details of this container are important. Daya and Ecke  find that the bulk temperature and velocity fluctuations in a cylindrical cell and in a square cell are very different from each other. For cylindrical cells the aspect ratio dependence has been examined by the Ahlers group  and by the Tong group ; in those experiments only the prefactors of effective power laws seem to change. Niemela and Sreenivasan  find an aspect ratio dependence of the heat flow at large Rayleigh numbers, and conclude that the focus of the next generation of experiments must be on large aspect ratio cells.

Another recent surprising result to be mentioned in this context is the dependence of the kinetic boundary layer (BL) thickness (which is set by the mean flow velocity profile above the wall) on its location in the cell. Qiu and Xia  find that it scales differently at the sidewalls and at the top or bottom plates. They find $`\delta _w/L=3.6Ra^{0.26\pm 0.03}`$ for the sidewall BL thickness (from now on abbreviated with w, for wall) and $`\delta _p/L=0.65Ra^{0.16\pm 0.02}Pr^{0.24\pm 0.01}`$ for the thickness at the top or bottom plates (from now on abbreviated with p, for plate), both for an aspect ratio $`\mathrm{\Gamma }=1`$ cylindrical cell. No existing theory of Rayleigh-Bénard convection accounts for this difference. Our own unifying theory of thermal convection  so far allows only for one kinetic boundary layer thickness, and gives a power law exponent $`0.22`$ (for $`\mathrm{\Gamma }=1`$ and a Prandtl number of $`Pr=6`$) , close to the average $`(0.160.26)/2=0.21`$ of the experimental results for the plate and the wall BLs thicknesses.

In this paper we set out to offer a mechanism to account for the shape and aspect ratio dependence in Rayleigh-Bénard convection, and for the difference between the wall and plate kinetic BLs. The starting point is Prandtl’s BL theory . The main physical ingredient is the conservation of the volume flux of the large scale wind. In the second part of the paper we present scaling arguments for the fluctuations of the velocity and the temperature in the bulk.

Laminar BL flow is described by the famous Prandtl equations . In these equations streamwise lengths are scaled by the streamwise length scale $`l`$, while wall-normal lengths are (re)scaled by $`l/\sqrt{Re}`$. Here, $`Re=lU_0/\nu `$ is the Reynolds number based on the streamwise length scale $`l`$ and the streamwise velocity scale $`U_0`$. Correspondingly, streamwise velocities are scaled by $`U_0`$ and wall-normal velocities by $`U_0/\sqrt{Re}`$. The immediate consequence is that the thickness of the laminar BL scales as $`\delta l/\sqrt{Re}`$; see also e.g. section 39 of ref. ).

In the context of Rayleigh-Bénard convection this means that the relevant streamwise length scales are the width $`d`$ of the cell for the plates and the height $`L`$ of the cell for the walls, i.e., different for cells with non-unity aspect ratio $`\mathrm{\Gamma }=d/L`$. Correspondingly, according to the Prandtl theory, also the widths of the BLs at the walls and at the plates are different, namely,

$$\delta _w=a\frac{L}{\sqrt{U_wL/\nu }}$$
(1)
and

$$\delta _p=a\frac{d}{\sqrt{U_pd/\nu }},$$
(2)
respectively. Here, $`a`$ is a dimensionless prefactor of order 1; Blasius  gave $`a=1.72`$ for a semi-infinite kinetic BL, cf. also . If we assumed that the streamwise wall and plate velocities are the same, $`U_w=U_p`$, we would immediately get aspect ratio dependent BL widths $`\delta _w=\delta _p/\sqrt{\mathrm{\Gamma }}`$.

However, the experiments by Tong’s group  show that the assumption $`U_w=U_p`$ does not hold. Indeed, a more realistic assumption seems to us volume flux conservation, i.e., equal volume fluxes of the wind over the walls and the plates,

$$\dot{V}_w=\dot{V}_p.$$
(3)
Then the geometry of the cell immediately enters. We will discuss two extreme situations: The “confined” flow, and the “unconfined” flow, see figure 1. In the confined flow the convection roll has the same spanwise extension both at the plates and the walls. In contrast, the unconfined flow makes use of the full plate width $`d`$ in the middle, whereas towards the walls its extension is of order $`2\sqrt{d\delta _w}`$ . The real flow will presumably be in between these two extremes.

Confined flow: The spanwise length scale of the flow is of order $`2\sqrt{d\delta _w}`$ at the walls and at the plates (see figure 1). The corresponding volume fluxes are $`\dot{V}_w=2\sqrt{d\delta _w}\delta _wU_w`$ at the walls and $`\dot{V}_p=2\sqrt{d\delta _w}\delta _pU_p`$ at the plates. With volume flux conservation eq. (3) and the Prandtl relations (1) and (2) we immediately get

$$\delta _w/\delta _p=\mathrm{\Gamma }^1,$$
(4)
$$U_p/U_w=\mathrm{\Gamma }^1.$$
(5)
Thus the wall and plate BLs have different thicknesses. But there is no difference in the Ra-scaling of $`\delta _w`$ and $`\delta _p`$, in contrast to what was found by Xia’s group .

Unconfined flow: The spanwise length scale of the flow at the cylindrical walls still is $`2\sqrt{d\delta _w}`$. But now by definition the flow at the plates is unconfined spanwise and the respective spanwise extension is of order $`d`$. The corresponding volume fluxes are now $`\dot{V}_w=2\sqrt{d\delta _w}\delta _wU_w`$ at the wall and $`\dot{V}_p=d\delta _pU_p`$ at the plate. With volume flux conservation eq. (3) and the Prandtl relations (1) and (2) we now obtain $`U_p^2/U_w=\nu L^1\mathrm{\Gamma }^42^4a^4`$ and

$$\delta _w/\delta _p=\mathrm{\Gamma }^{1/2}\sqrt{U_p/U_w}.$$
(6)
We now introduce wall and plate Reynolds numbers,

$$Re_w=U_wL/\nu ,Re_p=U_pL/\nu .$$
(7)
Note that both are defined with the height $`L`$ of the cell. With these definitions we find

$$\delta _w/\delta _p=\mathrm{\Gamma }^{5/2}4aRe_p^{1/2}=\sqrt{4a}\mathrm{\Gamma }^{3/2}Re_w^{1/4},$$
(8)
$$U_p/U_w=Re_p/Re_w=\mathrm{\Gamma }^416a^2Re_p^1=4a\mathrm{\Gamma }^2Re_w^{1/2}.$$
(9)
Equation (8) is our first main result. It shows, for unconfined flow, indeed a Reynolds number and thus a Rayleigh number dependence of the ratio between the BL thicknesses at the wall and at the plate, due to the cylindrical shape of the container. $`\delta _w`$ and $`\delta _p`$ scale differently with $`Ra`$, just as experimentally observed .

Since experimental flow is expected to lie in between the confined and the unconfined cases, $`\delta _w/\delta _p`$ should be between $`Re_p^{1/2}`$ and $`Re_p^0`$ or between $`Ra^{0.22}`$ and $`Ra^0`$; here we have used the effective scaling $`Re_pRa^{0.45}`$ which holds for $`Pr=5`$ and $`\mathrm{\Gamma }=1`$ both experimentally  and theoretically  in the $`Ra`$ regime $`10^8<Ra<10^{10}`$. Indeed, the experimental ratio has a scaling exponent in between the values given by the unconfined and confined cases, namely $`\delta _w/\delta _pRa^{0.11}`$ .

We now come to our other result, eq.(9), and to the aspect ratio dependence of $`\delta _w/\delta _p`$ and $`U_p/U_w`$. For $`\mathrm{\Gamma }=1/2`$ and $`\mathrm{\Gamma }=2`$ these ratios (as predicted in the two cases, and gauged to $`\mathrm{\Gamma }=1`$) are summarized in table I. The experimental data are obtained from figure 3 ($`\mathrm{\Gamma }=1`$, $`Ra=3.710^9`$), figure 6 ($`\mathrm{\Gamma }=2`$, $`Ra=4.910^8`$), and figure 8 ($`\mathrm{\Gamma }=1/2`$, $`Ra=3.2810^{10}`$) of reference . In this series of figures both $`\mathrm{\Gamma }`$ and $`Ra`$ were changed at the same time. In order to compare the ratios at $`\mathrm{\Gamma }=2`$ measured at $`Ra=510^8`$ with those at $`\mathrm{\Gamma }=1`$ measured at $`Ra=3.710^9`$, they first have to be extrapolated to this higher $`Ra`$ number. Correspondingly, in the $`\mathrm{\Gamma }=1/2`$ case we have to extrapolate the ratios from the original value $`Ra=3.310^{10}`$ to $`Ra=3.710^9`$. These extrapolations have been done according to the unconfined case formulas (8), (9) or according to the confined case ones (4), (5). In the latter case in fact no extrapolation is necessary as in those equations there is no $`Ra`$ or $`Re`$ dependence.

From table I we conclude that the experimental results (if extrapolated with the unconfined flow model) in general lie in between the predictions of the confined and unconfined cases, just as expected. An exception is $`U_p/U_w`$ for $`\mathrm{\Gamma }=2`$; we have no explanation for that. If extrapolated with the confined flow model both velocity ratios do not lie in between the predictions of the two flow cases.

The physics of the aspect ratio dependence is as follows: For increasing aspect ratio the available plate area of the flow increases and the plate velocity can go down. According to the Prandtl law (2), the plate BL thickness then increases. At the walls there is no additional lateral space for larger aspect ratios, therefore the flow has to accelerate, $`U_w>U_p`$, and correspondingly $`\delta _w<\delta _p`$.

Note that above dependence on $`\mathrm{\Gamma }`$ was derived assuming there is one convection roll. Once the convection develops two rolls on top of each other (for $`\mathrm{\Gamma }`$ less than 1) or next to each other (for $`\mathrm{\Gamma }`$ larger than 1), the aspect ratio dependence has to be embodied in a more sophisticated way. For small $`\mathrm{\Gamma }1/2`$ a transition from one roll to two rolls has been observed numerically . Possibly also in experiment both states of convection have been realized, see the discussion in . Thus one cannot formally consider $`\mathrm{\Gamma }0`$ or $`\mathrm{\Gamma }\mathrm{}`$ in the relations (4) – (9).

Can one manipulate a system such that it is no longer in between the confined and the unconfined cases? An option may be to take a square box and slightly tilt it in order to disfavor the flow along the diagonal, which would occur without tilting . If this can be achieved, one expects such flow to follow the confined flow model with confinement by the left and right walls. For this type of flow eq. (4) should hold, i.e., the wall and the plate BLs should show the same $`Ra`$-scaling. An experimental test of this predictions seems worthwhile to us to check the proposed theory.

Kinetic and thermal dissipation in the BLs: Within above framework the aspect ratio dependence can also be embodied into the unifying theory of thermal convection . The main idea of that theory is to split the kinetic dissipation $`ϵ_u`$ and the thermal dissipation $`ϵ_\theta `$, for which exact global relations can be derived from the Boussinesq equations, into their bulk and BL contributions.

For the kinetic dissipation this splitting obviously has to be extended to

$$ϵ_u=ϵ_{u,p}+ϵ_{u,w}+ϵ_{u,b},$$
(10)
where $`ϵ_{u,p}`$, $`ϵ_{u,w}`$ mean the kinetic dissipation in the plates’ or walls’ BLs, and $`ϵ_{u,b}`$ the bulk kinetic dissipation. Following , the BL dissipations are estimated as

$$ϵ_{u,p}=\nu \frac{U_p^2}{\delta _p^2}\frac{\delta _p}{L}=\frac{\nu ^3}{L^4}\frac{1}{a\mathrm{\Gamma }^{1/2}}Re_p^{5/2}=\frac{\nu ^3}{L^4}2^5a^{3/2}\mathrm{\Gamma }^{9/2}Re_w^{5/4}$$
(11)
and

$$ϵ_{u,w}=\nu \frac{U_w^2}{\delta _w^2}\frac{\delta _w}{l}=\frac{\nu ^3}{L^4}\frac{\mathrm{\Gamma }^9}{2^{10}a^6}Re_p^5=\frac{\nu ^3}{L^4}\frac{1}{a\mathrm{\Gamma }}Re_w^{5/2}.$$
(12)
The bulk dissipation will be dealt with below.

For the thermal dissipation the partition into BLs and bulk is different. Though at least in ideal RB convection the sidewalls are perfectly isolated (no flux condition), thermal BLs can develop at the sidewalls, as observed both in numerical simulations and in experiment . However, within our unifying theory ref.  it is not necessary to distinguish between sidewall thermal dissipation and bulk thermal dissipation. Therefore, here we introduce the notation $`ϵ_{\theta ,p}`$ for the thermal dissipation in the plate BLs and $`ϵ_{\theta ,\overline{p}}`$ for the thermal dissipation elsewhere, i.e., within the bulk and within the sidewall BLs (“non-plate”).

To estimate these contributions we remind that $`ϵ_\theta =\kappa \mathrm{\Delta }^2L^2Nu`$ and that the heat current $`Nu`$ consists of two terms

$`Nu`$ $`=`$ $`\left(u_z\theta _{A,t}(z)\kappa _z\theta _{A,t}(z)\right)/(\kappa \mathrm{\Delta }L^1)`$ (13)
$`=`$ $`\left(ϵ_{\theta ,\overline{p}}+ϵ_{\theta ,p}\right)/(\kappa \mathrm{\Delta }^2L^2)=ϵ_\theta /(\kappa \mathrm{\Delta }^2L^2).`$ (14)
Here $`\mathrm{}_{A,t}`$ denotes the average on time and $`x,y`$-plane.

For $`z`$ outside of the plates’ thermal BLs the first term in eq. (14) dominates. We estimate $`u_z\theta _{A,t}U\mathrm{\Delta }`$ and obtain $`ϵ_{\theta ,\overline{p}}/(\kappa \mathrm{\Delta }^2/L^2)PrRe`$. Note again that the thermal dissipation in the sidewall BLs is included in $`ϵ_{\theta ,\overline{p}}`$ and may even be dominant as compared to the contributions from the bulk, as the numerical simulations by Verzicco and Camussi suggest .

For $`z`$ within the plates’ thermal BLs, i.e., $`z/L0`$ or $`z/L1`$, the second term in eq. (14) dominates. Estimating the temperature gradient as $`\mathrm{\Delta }/\lambda _\theta `$ gives $`ϵ_{\theta ,p}/(\kappa \mathrm{\Delta }^2/L^2)L/\lambda _\theta `$. To connect the BL width $`\lambda _\theta `$ with $`Re`$, we use the temperature equation $`_t\theta =u_j_j\theta +\kappa _j_j\theta `$, leading to $`Ud^1\kappa \lambda _\theta ^2`$, which implies $`ϵ_{\theta ,p}/(\kappa \mathrm{\Delta }^2/L^2)\sqrt{RePr}`$, cf. . According to this derivation the Reynolds number one has to use in $`ϵ_{\theta ,\overline{p}}`$ is $`Re_w`$, while $`ϵ_{\theta ,p}`$ is determined by $`Re_p`$.

Scaling relations for bulk velocity and temperature fluctuations: Up to now we dealt with the BLs. As for the bulk fluctuations $`u^{}`$, the experiments  suggest that they scale differently from the large scale velocity $`U`$, namely with a weaker $`Ra`$ dependence. Typically, the fluctuation Reynolds number is $`Re^{}=u^{}L/\nu Ra^{0.40\pm 0.03}`$ and the large scale velocity Reynolds number with $`Re=UL/\nu Ra^{0.43\mathrm{}0.50}`$ . In refs.  we developed a theory for the dependences $`Re(Ra,Pr)`$ and $`Nu(Ra,Pr)`$. Here we add how the fluctuations $`Re^{}`$ and $`\theta ^{}`$ should behave as functions of $`Ra`$ and $`Pr`$.

The turbulence in the bulk is driven by the large scale wind, therefore $`ϵ_{u,b}U^3/L`$ as used in refs. . The energy cascades down and will be dissipated on scales comparable to the Kolmogorov scale $`\eta =(\nu ^3/ϵ_{u,b})^{1/4}`$. Thus we can estimate the bulk dissipation also with $`ϵ_{u,b}\nu u^2/\eta ^2`$. Combining both we obtain

$$Re^{}Re^{3/4}.$$
(15)
With $`ReRa^{0.43\mathrm{}0.50}`$ this means $`Re^{}Ra^{0.32\mathrm{}0.38}`$, which is close to above quoted experimental finding $`Re^{}Ra^{0.40\pm 0.03}`$ by Lam et al.  and also to Daya and Ecke’s finding  $`Re^{}Ra^{0.36\pm 0.05}`$ for a square cell.

The same reasoning can be followed to calculate the scaling of the typical bulk temperature fluctuations $`\theta ^{}`$. The thermal bulk dissipation is driven by the large scale temperature $`\mathrm{\Delta }`$; its time scale again is $`L/U`$. Thus $`ϵ_{\theta ,b}\mathrm{\Delta }^2U/L`$. On the other hand, the thermal energy in the bulk is dissipated by the temperature fluctuations $`\theta ^{}`$ on scale $`\eta _\theta =(\kappa ^3/ϵ_{u,b})^{1/4}LPr^{3/4}Re^{3/4}`$, i.e., $`ϵ_{\theta ,b}\kappa \theta ^2/\eta _\theta ^2`$. Balancing these two expressions for $`ϵ_{\theta ,b}`$ one obtains

$$\theta ^{}/\mathrm{\Delta }Pr^{1/4}Re^{1/4}.$$
(16)
With the experimental result  $`ReRa^{0.43}Pr^{0.76}`$, which holds in $`10^8<Ra<310^{10}`$ and $`3<Pr<1200`$, this means $`\theta ^{}/\mathrm{\Delta }Pr^{0.06}Ra^{0.11}`$. For $`ReRa^{1/2}`$  one gets $`\theta ^{}/\mathrm{\Delta }Ra^{0.13}`$. Both power laws are very close to the experimental findings $`\theta ^{}/\mathrm{\Delta }Ra^{1/7}`$  or $`\theta ^{}/\mathrm{\Delta }Ra^{0.10\pm 0.02}`$  for cylindrical cells.

In conclusion, we have offered a possible mechanism to cope with the recently observed effects of the RB cell geometry and the different $`Ra`$ behaviors of the kinetic BLs near the plates and the sidewalls. This has to be confirmed in more experimental detail. In particular, the comparison of RB convection in a square cell with diagonal (unconfined case) and with edge parallel (confined case) wind seems promising. – Another aspect is to distinguish the 1-roll from the 2-rolls (or more) states. It seems also clear that measuring the large scale wind velocity should take notice that possibly $`U_pU_w`$. Circulation time measurements average over both, while local methods can distinguish $`U_p`$ and $`U_w`$. – Finally, we have indicated how the global theory of $`Nu`$ versus $`Ra,Pr`$ has to be modified to take care of geometry effects and how the different scalings of bulk fluctuations and large scale quantities are related.

Acknowledgment: We thank X. van Doornum for drawing figure 1. The work is part of the research program of FOM, which is financially supported by NWO. It was also supported by the European Union (EU) under contract HPRN-CT-2000-00162 and by the German-Israeli Foundation (GIF).


text:
# Formation of Binary Pulsars in Globular Clusters

## 1. Introduction

The properties of globular cluster pulsars (see articles by Lorimer et al. and Ransom et al. in this volume) are rather surprising. While some pulsars are single, the majority are in short-period binaries. Most of the binaries have properties similar to those of the rare “eclipsing binary pulsars” seen in the Galactic disk population (see Nice et al. 2000 for a review). These systems have extremely short orbital periods, $`P_b110`$hr, circular orbits, and very low-mass companions, with $`m_20.03M_{}`$. The other, “normal” binaries have properties more similar to those of the disk population of binary millisecond pulsars, with nearly-circular orbits, periods $`P_b1`$d (near the short-period end of the distribution for binaries in the disk) and white-dwarf (WD) companions with $`m_2\mathrm{sin}i0.2M_{}`$.

The large inferred total population of millisecond pulsars in globular clusters (several hundred in 47 Tuc alone; see Camilo et al. 2000) and the very high stellar densities in many cluster cores ($`\rho _c10^410^6M_{}\mathrm{pc}^3`$) suggest that dynamical interactions must play a dominant role in the formation of these systems. However, the dynamical formation scenarios traditionally invoked for the production of recycled pulsars in globular clusters have many problems.

Scenarios based on tidal capture of low-mass main-sequence (MS) stars by neutron stars (NS), followed by accretion and recycling of the NS during a stable mass-transfer phase, run into many difficulties. Serious problems have been pointed out about the tidal capture process itself (which, because of strong nonlinearities in the regime relevant to globular clusters, is far more likely to result in a merger than in the formation of a detached binary; see, e.g., Kumar & Goodman 1996; McMillan et al. 1990; Rasio & Shapiro 1991). Moreover, the basic predictions of tidal capture scenarios are at odds with many observations of binaries and pulsars in clusters (Bailyn 1995; Johnston et al. 1992; Shara et al. 1996). It is likely that “tidal-capture binaries” are either never formed, or contribute negligibly to the production of recycled pulsars. Physical collisions between NS and red giants have also been invoked as a way of producing directly NS–WD binaries with ultra-short periods (e.g., Verbunt 1987), but detailed hydrodynamic simulations show that this does not occur (Rasio & Shapiro 1991).

The viability of tidal capture and two-body collision scenarios has become less relevant with the realization over the last decade that globular clusters contain dynamically significant populations of primordial binaries (Hut et al. 1992). Dynamical interactions involving hard primordial binaries are now thought to provide the dominant energy production mechanism that allows many globular clusters to remain in thermal equilibrium and avoid core collapse over $`10^{10}`$yr (Gao et al. 1991; McMillan & Hut 1994; Fregeau et al. 2003).

In dense clusters, retained NS can acquire binary companions through exchange interactions with these primordial binaries. Because of its large cross section, this process dominates over any kind of two-body interaction even for low primordial binary fractions (Heggie, Hut, & McMillan 1996; Leonard 1989; Sigurdsson & Phinney 1993). In contrast to tidal capture, exchange interactions with hard primordial binaries (with semimajor axes $`a0.11`$AU) can form naturally the wide binary millisecond pulsars seen in some low-density globular clusters (such as PSR B1310$`+`$18, with $`P_b=256`$d, in M53, which has the lowest central density, $`\rho _c10^3M_{}\mathrm{pc}^3`$, of any globular cluster with observed radio pulsars). When the newly acquired MS companion, of mass $`<1M_{}`$, evolves up the giant branch, the orbit circularizes and a period of stable mass transfer begins, during which the NS is recycled (see, e.g., Rappaport et al. 1995). The resulting NS–WD binaries have orbital periods in the range $`P_b110^3`$d. However, this scenario cannot explain the formation of recycled pulsars in binaries with periods shorter than $`1`$d. To obtain such short periods, the initial primordial binary must be extremely hard, with $`a<0.01`$AU, but then the recoil velocity of the system following the exchange interaction would almost certainly exceed the escape speed from the shallow cluster potential (e.g., $`v_e60\mathrm{km}\mathrm{s}^1`$ for 47 Tuc).

One can get around this problem by considering more carefully the stability of mass transfer in NS–MS binaries formed through exchange interactions. While all MS stars in the cluster today have masses $`<1M_{}`$, the rate of exchange interactions should have peaked at a time when significantly more massive MS stars were still present. Indeed, the NS and the most massive primordial binaries will undergo mass segregation and concentrate in the cluster core on a time scale comparable to the initial half-mass relaxation time $`t_{\mathrm{rh}}`$. For typical dense globular clusters, we expect $`t_{\mathrm{rh}}10^9`$yr, which is comparable to the MS lifetime of a $`23M_{}`$ star. If the majority of NS in the cluster core acquire MS companions in the range of $`13M_{}`$, a very different evolution ensues. Indeed, in this case, when the MS star evolves and fills its Roche lobe, the mass transfer for many systems (depending on the mass ratio and evolutionary state of the donor star) is dynamically unstable and leads to a common-envelope (CE) phase (see, e.g., Taam & Sandquist 2000). The emerging binary will have a low-mass WD in a short-period, circular orbit around the NS.

This simple idea has been explored quantitatively using Monte Carlo simulations of globular cluster dynamics (Rasio, Pfahl, & Rappaport 2000; Rappaport et al. 2001) and a brief summary of this work is presented in Sec. 2. A similar scenario, but starting from tidal capture binaries and applied to X-ray sources in globular clusters, was discussed by Bailyn & Grindlay (1987). The possibility of forming intermediate-mass binaries through exchange interactions was mentioned by Davies & Hansen (1998), who pointed out that NS retention in globular clusters may also require that the NS be born in massive binaries. Among eclipsing pulsars in the disk, at least one system (PSR J2050$``$0827) is likely to have had an intermediate-mass binary progenitor, given its very low transverse velocity (Stappers et al. 1998).

## 2. Results from Dynamical Simulations

Figure 1 shows typical results from our most recent numerical simulations for 47 Tuc (Rappaport et al. 2001). Given the many simplifying assumptions used in this work, the agreement with observations is quite encouraging. The parameters of observed systems were taken from Camilo et al. (2000) and Freire et al. (2001, 2003). We omitted 47 Tuc W, which is now known to have a MS companion (Edmonds et al. 2002). Such systems are the natural descendants of longer-period binaries (upper part of group c in Fig. 1), which are very likely to interact again in the dense cluster environment. These multiple interactions can produce both single recycled pulsars and binary pulsars with anomalous companions (like 47 Tuc W) and/or unusual positions in the cluster halo (as in NGC 6752; D’Amico et al. 2002).

The initial conditions for the simulations include MS stars and primordial binaries with standard initial mass functions and distributions of binary parameters, as well as a certain number of single NS, which are assumed to have been retained by the cluster (cf. Pfahl, Rappaport, & Podsiadlowski 2002). The total number of MS stars initially is $`10^610^7`$ and the total number of NS is $`10^4`$.

Binaries and NS undergo mass segregation and enter the cluster core in a time $`t_\mathrm{s}10(m_f/m_t)t_{\mathrm{rh}}`$, for objects of mass $`m_t`$ drifting through field stars of average mass $`m_f`$ (see, e.g., Fregeau et al. 2000). For the simulation of Fig. 1 we assumed $`t_{\mathrm{rh}}=10^9`$yr. Binaries whose primaries evolve off the MS before entering the cluster core are removed from the simulation. From the numbers of binaries and NS in the core, we can compute the time for each binary to have a strong interaction with a NS. We can then decide on a list of actual interactions in each timestep using a Monte Carlo procedure. All 3–body interactions are computed by direct numerical integration. Disrupted and ejected binaries are removed from the simulation.

We then calculate the evolution of the newly formed NS–MS binaries. When the primary evolves off the MS, the orbit is assumed to circularize (conserving total angular momentum). We then test for the stability of mass transfer when the primary fills its Roche lobe. We find that typically about 50% of the systems enter a CE phase. The outcome of the CE phase is calculated using the standard treatment, with the efficiency parameter $`\alpha _{CE}=0.5`$ (Rappaport, Di Stefano, & Smith 1994). A significant fraction of these NS–WD binaries will undergo further evolution driven by gravitational radiation. For orbital periods $`<8`$hr, the companion will be filling its Roche lobe in less than $`10^{10}`$yr and a second phase of mass transfer will occur. For WD masses $`<0.4M_{}`$ the mass transfer is stable and the evolution can be calculated semi-analytically. Our calculations for this phase also incorporate a simple treatment of the tidal heating of the companion (Applegate & Shaham 1994; Rasio et al. 2000). We track the accretion rate and spin-up of the NS during the mass-transfer phase and we terminate the evolution when the NS spin period reaches a randomly chosen value in the range $`25`$ms (at which point the radio pulsar emission is assumed to turn on and stop the accretion flow).

This simple scenario provides a natural way of explaining the large number and observed properties of short-period binary pulsars in a dense globular cluster such as 47 Tuc (Fig. 1). Although quantitatively the predicted properties of the final binary population depend on our parametrization of several uncertain processes (such as CE evolution and tidal heating), the overall qualitative picture is remarkably robust. Indeed, quite independent of the details of the various assumptions and choices of parameters, exchange interactions inevitably form a large population of NS–MS binaries that will go through a CE phase. The only way for a globular cluster to avoid forming such a population would be to start with a very low primordial binary fraction, a very small number of retained NS, or to have a very long relaxation time $`t_{\mathrm{rh}}>10^{10}`$yr, such that all MS stars with masses $`>1M_{}`$ evolve before the rate of exchange interactions becomes significant. A large fraction of the post-CE NS–WD binaries cannot avoid further evolution driven by gravitational wave emission, with the companion ultimately reduced to a very low mass $`m_210^2M_{}`$.

Pulsars in wider binaries with $`P_b1`$d and WD companions must have evolved from the group of systems with stable mass transfer from a $`1M_{}`$ subgiant to a NS, which have orbital periods at the start of mass transfer in the range $`15`$d (lower end of group c in Fig. 1). Conventional evolutionary scenarios suggest that systems where the donor has a well-developed degenerate core should inevitably evolve to longer orbital periods. However, many of the systems in group c of Fig. 1 have not yet developed such cores, and their final fate is uncertain. We also note that many binary pulsars in the Galactic disk population, which are supposed to fit this evolutionary scenario involving stable mass transfer from a low-mass giant to the NS, have orbital periods shorter than $`5`$d, with some $`<1`$d. Even the most sophisticated binary evolution models currently available cannot explain these systems (see Podsiadlowski, in this volume).

The results illustrated in Fig. 1 also predict the existence of a large number of binary pulsars with companion masses $`m_20.030.05M_{}`$ and orbital periods as short as $`15`$min that may have so far escaped detection (lower end of group a in Fig. 1). Future observations using more sophisticated acceleration-search techniques or shorter integration times may be able to detect them.

### Acknowledgments.

I thank Eric Pfahl and Saul Rappaport for many important contributions to this work. This work was supported by NASA Grants NAG5-8460, NAG5-11396, and NAG5-12044.

## References

Applegate, J.H., & Shaham, J. 1994, ApJ, 436, 312

Bailyn, C.D. 1995, ARAA, 33, 133

Camilo, F., Lorimer, D.R., Freire, P., Lyne, A.G., & Manchester, R.N. 2000, ApJ, 535, 975

D’Amico, N., et al. 2002, ApJ, 570, L89

Davies, M.B., & Hansen, B.M.S. 1998, MNRAS, 301, 15

Edmonds, P.D., Gilliland, R.L., Camilo, F., Heinke, C.O., Grindlay, J.E. 2002, ApJ, 579, 741

Fregeau, J., Joshi, K., Portegies Zwart, S., & Rasio, F.A. 2002, ApJ, 570, 171

Fregeau, J.M., Gürkan, A., Joshi, K.J., & Rasio, F.A. 2003, ApJ, submitted

Freire, P.C., et al. 2001, MNRAS, 326, 901

Freire, P.C., et al. 2003, MNRAS, submitted

Gao, B., Goodman, J., Cohn, H., & Murphy, B. 1991, ApJ, 370, 567

Johnston, H.M., Kulkarni, S.R., & Phinney, E.S. 1992, in X-Ray Binaries and Recycled Pulsars, eds. E.P.J. van den Heuvel & S.A. Rappaport (Dordrecht: Kluwer), 349

Heggie, D.C., Hut, P., & McMillan, S.L.W. 1996, ApJ, 467, 359

Hut, P., et al. 1992, PASP, 104, 981

Kumar, P., & Goodman, J. 1996, ApJ, 466, 946

Leonard, P.J.T. 1989, AJ, 98, 217

McMillan, S., & Hut, P. 1994, ApJ, 427, 793

McMillan, S.L.W., Taam, R.E., & McDermott, P.N. 1990, ApJ, 354, 190

Nice, D.J., Arzoumanian, Z., & Thorsett, S.E. 2000, in IAU Colloq. 177, Pulsar Astronomy - 2000 and Beyond, eds. M. Kramer et al. (ASP Conf. Series Vol. 202), 67

Pfahl, E., Rappaport, S., & Podsiadlowski, P. 2002, ApJ, 573, 282

Rappaport, S., Di Stefano, R., & Smith 1994, ApJ, 426, 692

Rappaport, S., Pfahl, E., Rasio, F.A., & Podsiadlowski, P. 2001, in Evolution of Binary and Multiple Star Systems, eds. P. Podsiadlowski et al. ASP Conf. Series, Vol. 229), 409

Rappaport, S., Podsiadlowski, P., Joss, P.C., Di Stefano, R., & Han, Z. 1995, MNRAS, 273, 731

Rasio, F.A., Pfahl, E.D., & Rappaport, S.A. 2000, ApJ, 532, L47

Rasio, F.A., & Shapiro, S.L. 1991, ApJ, 377, 559

Shara, M.M., Bergeron, L.E., Gilliland, R.L., Saha, A., & Petro, L. 1996, ApJ, 471, 804

Sigurdsson, S., & Phinney, E.S. 1993, ApJ, 415, 631

Stappers, B.W., Bailes, M., Manchester, R.N., Sandhu, J.S., & Toscano, M. 1998, ApJ, 499, L183

Taam, R.E., & Sandquist, E.L. 2000, ARAA, 38, 113

Verbunt, F. 1987, ApJ, 312, L23


text:
# Some remarks on cyclization and periodicity

## 1 Introduction

Let us put in correspondence to any function $`F\left(\phi \right)`$ function

$$\stackrel{~}{F}\left(\phi \right)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}F\left(\phi +2\pi n\right)$$
(1)
This procedure $`\left(\text{1}\right)`$, which we call cyclization, is used not infrequently (see e.g. ), but as a rule is applied to such functions $`F\left(\phi \right)`$ that decrease quite rapidly with $`\phi \pm \mathrm{}`$ and $`\stackrel{~}{F}\left(\phi \right)`$ turns out to be periodic. In this paper we consider cases when the above condition is not true and $`\stackrel{~}{F}\left(\phi \right)\stackrel{~}{F}\left(\phi +2\pi \right)`$.

It is worth noting from the outset that since shift $`\stackrel{~}{F}\left(\phi \right)`$ $`\stackrel{~}{F}\left(\phi +2\pi \right)`$ leads to some sort of interchange of summations, function $`\stackrel{~}{F}\left(\phi \right)`$ may appear to be aperiodic, if series $`\left(\text{1}\right)`$ is not absolutely convergent.

Indeed, if we take as an example

$$F\left(\phi \right)=\mathrm{ln}\left(1+\frac{\alpha }{\phi }\right)\frac{\alpha }{\phi },$$
(2)
then after cyclization we get

$$\stackrel{~}{F}\left(\phi \right)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{ln}\left(\left(1+\frac{\alpha }{\phi +2\pi n}\right)\mathrm{exp}\left\{\frac{\alpha }{\phi +2\pi n}\right\}\right)$$
(3)
which, taking into account ( 6.2.3.$`\left(8\right)`$) may be written as

$$\stackrel{~}{F}\left(\phi \right)=\frac{\alpha }{2\pi }\mathrm{Psi}\left(\frac{\phi }{2\pi }\right)+\mathrm{ln}\frac{\mathrm{\Gamma }\left(\frac{\phi }{2\pi }\right)}{\mathrm{\Gamma }\left(\frac{\phi +\alpha }{2\pi }\right)}$$
(4)
It is easy to check that

$$\stackrel{~}{F}\left(\phi +2\pi \right)\stackrel{~}{F}\left(\phi \right)=\frac{\alpha }{\phi }\mathrm{ln}\left(1+\frac{\alpha }{\phi }\right)0$$
(5)
In many interesting cases series in $`\left(\text{1}\right)`$ is divergent in ordinary sense, then to compute $`\stackrel{~}{F}\left(\phi \right)`$ it is rational to apply distribution theory where rigorous methods is developed for summing extensive class of such series . In particular, such methods allow to compute integrals of tempered distributions and series which terms coincide with some tempered function at integer values of the argument .

Recall that $`F\left(\phi \right)`$ is tempered distribution, if and only if it is finite order derivative of some tempered function. Continuous $`f\left(\phi \right)`$ function is called a tempered one , if exists some $`\sigma `$ that

$$\left|f\left(\phi \right)\right|<\left|\phi \right|^\sigma ;\left|\phi \right|\mathrm{}.$$
(6)
It is important that every periodic distribution $`\stackrel{~}{F}_{\left(p\right)}\left(\phi \right)`$ is tempered  and Fourier transform of any tempered distribution is tempered distribution as well .

It is known that any tempered distribution $`F\left(\phi \right)`$ or $`\stackrel{~}{F}_{\left(p\right)}\left(\phi \right)`$ may be presented as (see e.g. ).

$$F\left(\phi \right)=F_+\left(\phi +i\epsilon \right)F_{}\left(\phi i\epsilon \right)$$
(7)
and

$$\stackrel{~}{F}_{\left(p\right)}\left(\phi \right)=\stackrel{~}{F}_+\left(\phi +i\epsilon \right)\stackrel{~}{F}_{}\left(\phi i\epsilon \right);$$
(8)
where $`\stackrel{~}{F}_\pm \left(\phi \right)`$ and $`F_\pm \left(\phi \right)`$ are analytical functions in upper/lower complex half-plane $`\phi `$ and $`\epsilon `$ is routine positive infinitesimal parameter.

It may be easily checked, if to choose

$$F_+\left(\phi \right)_0^{\mathrm{}}F_te^{i\phi t}𝑑t;F_{}\left(\phi \right)_{\mathrm{}}^0F_te^{i\phi t}𝑑t$$
(9)
and

$$\stackrel{~}{F}_+\left(\phi \right)\underset{n=0}{\overset{\mathrm{}}{}}F_ne^{i\phi n};\stackrel{~}{F}_{}\left(\phi \right)\underset{n=\mathrm{}}{\overset{1}{}}F_ne^{i\phi n}$$
(10)
then expression $`\left(\text{7}\right)`$ will coincide with an Abel-Poisson regularization of Fourier integral

$$F\left(\phi \right)=_{\mathrm{}}^{\mathrm{}}F_te^{it\phi }𝑑t_{\mathrm{}}^{\mathrm{}}F_te^{it\phi \epsilon \left|t\right|}𝑑t$$
(11a)
and $`\left(\text{8}\right)`$ with an Abel-Poisson regularization of Fourier series
$$\stackrel{~}{F}_{\left(p\right)}\left(\phi \right)\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}F_ne^{i\phi n}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}F_ne^{i\phi n\epsilon \left|n\right|}$$
(12)
Abel-Poisson regularization provide continuous convergence Fourier integrals $`\left(\text{11a}\right)`$ and series $`\left(\text{12}\right)`$ in a case when $`F_t`$ is tempered distribution and Fourier coefficients $`F_n`$ coincide with values of some tempered distribution at integer $`t=n`$.

In particular, due to continuous convergence of regularized integrals, the order of integration in double integrals may be changed and it may be shown, that the transform inverse to $`\left(\text{11a}\right)`$ is

$$F_t=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}F\left(\phi \right)\mathrm{exp}\left\{i\phi t\epsilon \left|\phi \right|\right\}𝑑\phi ,$$
(13a)
Furthermore, since regularized111Here $`\epsilon ^{}`$ is positive infinitesimal parameter independent of $`\epsilon `$. sums and integrals in
$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}F\left(\phi 2\pi n\right)`$
$`={\displaystyle \underset{n=\mathrm{}}{\overset{1}{}}}{\displaystyle _{\mathrm{}}^0}e^{i\left(\phi 2\pi ni\epsilon \right)\left(t+i\epsilon ^{}\right)}F_t𝑑t+{\displaystyle \underset{n=\mathrm{}}{\overset{1}{}}}{\displaystyle _0^{\mathrm{}}}e^{i\left(\phi 2\pi n+i\epsilon \right)\left(t+i\epsilon ^{}\right)}F_t𝑑t+`$
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle _{\mathrm{}}^0}e^{i\left(\phi 2\pi ni\epsilon \right)\left(ti\epsilon ^{}\right)}F_t𝑑t+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}e^{i\left(\phi 2\pi n+i\epsilon \right)\left(ti\epsilon ^{}\right)}F_t𝑑t`$ (14)
are continuously convergent it allows to change the order of integration and summation, thus taking into account the Poisson relation

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{i2\pi nm}=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\delta \left(nm\right)$$
(15)
one can easily find that application of the procedure $`\left(\text{1}\right)`$ to regularized integral $`\left(\text{11a}\right)`$ gives

$$\stackrel{~}{F}\left(\phi \right)=_{\mathrm{}}^{\mathrm{}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{i\left(\phi 2\pi n\right)t}F_tdt=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}e^{i\phi t}\delta \left(tm\right)F_t𝑑t$$
(16)
In a case when for all integer $`n`$ function $`F_t`$ obeys the condition

$$F_n=\underset{tn}{lim}F_t<\mathrm{}$$
(17)
series $`\left(\text{16}\right)`$ coincides with $`\left(\text{12}\right)`$ and $`F_n`$ can be considered Fourier coefficients which may be computed, as it follows from $`\left(\text{1}\right)`$ and $`\left(\text{13a}\right)`$, from the ordinary relation

$$F_n=\frac{1}{2\pi }_\pi ^\pi \stackrel{~}{F}\left(\phi \right)\mathrm{exp}\left\{i\phi n\right\}𝑑\phi $$
(18)
For $`F_t`$ to exist function $`F\left(\phi \right)`$ should be only locally summable rather than summable. However, if $`F\left(\phi \right)`$ is summable, its Fourier transform $`F_t`$ is bounded for any real $`t`$ .

It is worth to note, that too rigid constraint on class function to which $`F\left(\phi \right)`$ belongs, may lead to quite trivial class function $`F_t`$ and, consequently, trivial $`\stackrel{~}{F}\left(\phi \right)`$. Indeed, let $`F\left(\phi \right)`$ be infinitely differentiable with compact support (i.e. belongs to Schwartz space $`\left(D\right)`$), then, according to , its Fourier transform $`F_t`$ belongs to space of entire functions $`\left(Z\right)`$. If in addition

$$\underset{\left|t\right|\mathrm{}}{lim}\frac{\mathrm{ln}\left|F_t\right|}{\left|t\right|}=0$$
(19)
for any direction in the complex plane $`t`$, and $`F_t`$ is bounded for all integer $`t=n`$, then, according Polya theorem (see e.g. ), $`F_n=const`$.

It is evident that, if condition $`\left(\text{17}\right)`$ is fulfilled, function $`\stackrel{~}{F}\left(\phi \right)`$ given by $`\left(\text{1}\right)`$ may be presented as regularized Fourier series $`\left(\text{12}\right)`$ and consequently is periodic.

## 2 Singular $`F_t`$

In a case when function $`F\left(\phi \right)`$ is not summable, Fourier transform $`F_t`$ may appear to be singular at integer $`t=n`$, as in already considered example $`\left(\text{2}\right)`$ where $`F_t`$, computed with (  402), may be written as

$$F_t=\left(e^{i\alpha t}1\right)\left(\pi \left(2\mathrm{Psi}\left(1\right)+1\right)\delta \left(t\right)\frac{i}{t}\frac{\pi }{\left|t\right|}\right)i\pi \alpha signum\left(t\right)$$
(20)
The case of singular $`F_t`$ needs special consideration. As it is seen from $`\left(\text{16}\right)`$, to obtain $`\stackrel{~}{F}\left(\phi \right)`$ we have to compute the product of two distributions: $`\delta \left(tn\right)`$ and $`F_t`$. Unfortunately, the product of arbitrary distributions is not defined (see e.g. ), so we shall confine ourself to the consideration of well-defined products, namely the case when singularity of $`F_t`$ at some integer $`t=k`$ is a pole of order $`m`$. More general case will be considered elsewhere. To be specific we consider $`k=0`$ i.e. $`F_tt^m`$ for $`t0`$.

Let us take as an example $`F\left(\phi \right)=signum\left(\phi \right)`$. Making allowance for ( 13), one may write

$$F\left(\phi \right)=signum\left(\phi \right)=\frac{1}{\pi i}_{\mathrm{}}^{\mathrm{}}e^{it\phi }t^1𝑑t$$
(21)
and in this case

$$F_t=\frac{1}{\pi i}t^1$$
(22)
is infinite at $`t=0`$.

It is easy to check that

$$\stackrel{~}{F}\left(\phi \right)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}signum\left(\phi +2\pi n\right)=\frac{1}{\pi i}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}\delta \left(mt\right)e^{it\phi }t^1𝑑t$$
(23)
which leads to

$$\stackrel{~}{F}\left(\phi \right)=\stackrel{~}{F}_{\left(p\right)}\left(\phi \right)+\stackrel{~}{F}_{\left(a\right)}\left(\phi \right)$$
(24)
where the first term is evidently periodic due to explicit invariance under $`\phi \phi +2\pi `$

$$\stackrel{~}{F}_{\left(p\right)}\left(\phi \right)=\underset{m0}{}\frac{e^{im\phi }}{m}=\frac{1}{\pi i}\mathrm{ln}\frac{1e^{i\phi }}{1e^{i\phi }}$$
(25)
Since a logarithm is a many-valued function, we get the ambiguous result

$$\frac{1}{\pi i}\mathrm{ln}\frac{1e^{i\phi }}{1e^{i\phi }}=\frac{1}{\pi i}\mathrm{ln}\left(e^{i\phi }\right)=\frac{1}{\pi }\left(\left(\phi \pm \pi \right)+2\pi n_\phi \right)$$
(26)
where $`n_\phi `$ is an integer number, which numbers the branches of the logarithm.

To reduce ambiguity, we choose, as it is generally done, among the infinite set of logarithm branches, the principal one

$$\frac{1}{\pi i}\underset{m0}{}\frac{e^{im\phi }}{m}=\frac{\left(\phi \pm \pi \right)_{\mathrm{mod}2\pi }}{\pi };$$
(27)
where the periodic function $`\left(x+2\pi \right)_{\mathrm{mod}2\pi }=x_{\mathrm{mod}2\pi }`$ is defined in the interval $`\pi <x<\pi `$ as $`x_{\mathrm{mod}2\pi }=x`$ .

The second term in $`\left(\text{24}\right)`$

$$\stackrel{~}{F}_{\left(a\right)}\left(\phi \right)=\frac{1}{\pi i}_{\mathrm{}}^{\mathrm{}}e^{it\phi }\delta \left(t\right)t^1𝑑t$$
(28)
has the singular integrand that with the help of a simple equation

$$\delta \left(t\right)t^nn!=\left(1\right)^n\delta ^{\left(n\right)}\left(t\right)$$
(29)
may be rewritten as

$$\stackrel{~}{F}_{\left(a\right)}\left(\phi \right)=\frac{1}{\pi i}_{\mathrm{}}^{\mathrm{}}e^{it\phi }\delta ^{\left(1\right)}\left(t\right)𝑑t=\frac{\phi }{\pi }$$
(30)
This term is evidently aperiodic.

Taking into account

$$\frac{\phi \phi _{\mathrm{mod}2\pi }}{2\pi }=\frac{\phi \pi }{2\pi }+1=\frac{\phi +\pi }{2\pi }$$
(31)
where $`x`$ stands for integer part of $`x`$ (called also $`Entier\left(x\right)`$), one may finally write

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}signum\left(\phi +2\pi n\right)=2\frac{\phi }{2\pi }+1$$
(32)
A more general case may be considered, making allowance ( 22)

$$\left(1\right)^m\frac{\left(2m\right)!}{\pi i}_{\mathrm{}}^{\mathrm{}}e^{it\phi }t^{2m1}𝑑t=\left|\phi \right|^{2m}signum\left(\phi \right)F\left(\phi \right)$$
(33)
that leads to

$`\stackrel{~}{F}\left(\phi \right)`$ $`={\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\left|\phi +2\pi n\right|^{2m}signum\left(\phi +2\pi n\right)`$ (34)
$`=\left(1\right)^m{\displaystyle \frac{\left(2m\right)!}{\pi i}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{it\left(\phi +2\pi n\right)}t^{2m1}𝑑t`$
so one may write

$$\stackrel{~}{F}\left(\phi \right)=\frac{\phi ^{2m+1}}{\left(2m+1\right)\pi }+\frac{\left(2\pi \right)^{2m+1}}{\pi }\left(2m\right)!\underset{k0}{}\frac{e^{ik\phi }}{\left(2\pi ik\right)^{2m+1}}$$
(35)
Recall that (see e.g.  1.13$`\left(11\right)`$)

$$n!\underset{k0}{}\frac{\mathrm{exp}\left\{ik\phi \right\}}{\left(2\pi ik\right)^n}=B_n\left(\frac{\phi }{2\pi }\right);$$
(36)
where $`B_n\left(x\right)`$ is Bernoulli polynomial of order $`n`$. Unfortunately, expression $`\left(\text{36}\right)`$ is true only for $`0<\phi <2\pi `$ . One may see, however, that an extension of $`\left(\text{36}\right)`$ on negative $`\phi `$ may be easily done by the substitution $`\phi \phi +2\pi `$, which doesn’t effect explicitly periodic right side of $`\left(\text{36}\right)`$, but transforms condition $`0<\phi <2\pi `$ into $`2\pi <\phi <0`$, and we may finally write

$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\left|\phi +2\pi n\right|^{2m}signum\left(\phi +2\pi n\right)`$ (37)
$`={\displaystyle \frac{\phi ^{2m+1}}{\left(2m+1\right)\pi }}{\displaystyle \frac{\left(2\pi \right)^{2m+1}}{\left(2m+1\right)\pi }}B_{2m+1}\left({\displaystyle \frac{\phi _{\mathrm{mod}2\pi }+2\pi \theta \left(\phi _{\mathrm{mod}2\pi }\right)}{2\pi }}\right).`$
Proceeding as in previous case and taking into account ( 21)

$$\left|\phi \right|^{2m1}=\left(1\right)^m\frac{\left(2m1\right)!}{\pi }_{\mathrm{}}^{\mathrm{}}e^{it\phi }\frac{1}{t^{2m}}𝑑t$$
(38)
we may compute

$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}\left|\phi +2\pi n\right|^{2m1}`$ $`=\left(1\right)^m{\displaystyle \frac{\left(2m1\right)!}{\pi }}{\displaystyle \underset{k0}{}}e^{ik\phi }k^{2m}`$ (39)
$`+\left(1\right)^m{\displaystyle \frac{\left(2m1\right)!}{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{it\phi }\delta \left(t\right)t^{2m}𝑑t`$
that leads to

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\left|\phi +2\pi n\right|^{2m1}=\frac{\left(2\pi \right)^{2m1}}{m}B_{2m}\left(\frac{\phi _{\mathrm{mod}2\pi }}{2\pi }+\theta \left(\phi \right)\right)+\frac{\phi ^{2m}}{2\pi m}$$

In particular for $`m=1`$ we get

$$\stackrel{~}{F}\left(\phi \right)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\left|\phi +2\pi n\right|=\frac{\phi ^2}{2\pi }\frac{\phi _{\mathrm{mod}2\pi }^2}{2\pi }+\left|\phi _{\mathrm{mod}2\pi }\right|\frac{\pi }{3}.$$

It should be noted that for some functions $`F\left(\phi \right)`$ after lifting Abel-Poisson regularization $`\stackrel{~}{F}\left(\phi \right)`$ ’converges’ to infinity. In particular, for $`F\left(\phi \right)=\phi ^me^{i\alpha \phi }`$

$$\stackrel{~}{F}\left(\phi \right)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\left(\phi +2\pi n\right)^me^{i\alpha \left(\phi +2\pi n\right)}=i^m\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{in\phi }\delta ^{\left(m\right)}\left(\alpha n\right)$$
(40)
We see that for noninteger $`\alpha `$ function $`\stackrel{~}{F}\left(\phi \right)=0`$ whereas for integer $`\alpha =k`$ one gets $`\stackrel{~}{F}\left(\phi \right)=\left(\phi \right)^me^{ik\phi }\delta \left(0\right)`$

## 3 Conclusions

This paper is an attempt to reveal a mechanism that actualizes a known connection between discretization of Fourier variables and periodicity in the associated space.

It is shown that the application of procedure $`\left(\text{1}\right)`$ to function $`F\left(\phi \right)`$ doesn’t lead to the periodicity of function $`\stackrel{~}{F}\left(\phi \right)`$ if Fourier transform of $`F\left(\phi \right)`$ has a pole at some integer value of its argument. If, on the contrary, $`F_t`$ is finite for all integer $`t=n`$, function $`\stackrel{~}{F}\left(\phi \right)`$ is periodic.


text:
# Two-dimensional projections of an hypercube

## Abstract

We present a method to project a hypercube of arbitrary dimension on the plane, in such a way as to preserve, as well as possible, the distribution of distances between vertices. The method relies on a Montecarlo optimization procedure that minimizes the squared difference between distances in the plane and in the hypercube, appropriately weighted. The plane projections provide a convenient way of visualization for dynamical processes taking place on the hypercube.

Dynamical models where the state of a system is represented by an ordered array of binary variables are ubiquitous in statistical physics, especially in its interdisciplinary applications. Perhaps the widest class of models that admit such Boolean-like representation is constituted by binary cellular automata cel ; kauf0 . Specific applications include biological evolution at the levels of molecules volkenstein , cells kauf , individuals Penna  and species KJ ; SM , as well as social and socioeconomical behavior Palmer ; Weis . Moreover, genetic algorithms are typically applied to systems whose configuration is described by means of binary sequences GA . These models may involve large populations of interacting agents, each of them described as a time-dependent array of bits, which requires assigning an evolving density to each possible binary sequence.

While the configuration space of a binary sequence of length $`L`$ is naturally represented as the set of $`2^L`$ vertices of an $`L`$-dimensional hypercube, its visualization can be dissapointingly difficult, even for $`L`$ not very large. On the other hand, besides a quantitative characterization of the system dynamics through its collective properties, it is sometimes desirable to rely on a geometrical depiction where the dynamics can be followed, for instance, on the computer screen. The purpose of this paper is to present a method to project the vertices of a hypercube of arbitrary dimension onto a set of points in the plane, with the condition of preserving, as much as possible, the structure of the distance distribution on the hypercube. The motivation of this condition is that many dynamical processes depend on the Hamming distance—i.e., the number of different bits—between binary sequences, and we require this feature to be well represented by the Euclidean distance between the corresponding points in the plane projection.

Let $`h_{ij}`$ be the Hamming distance between vertices $`i`$ and $`j`$ in the hypercube, and $`d_{ij}`$ the Euclidean distance between points $`i`$ and $`j`$ in their plane projection. We define the function

$$E=\underset{i,j}{}(d_{ij}h_{ij})^2,$$
(1)
that characterizes how different are the distances between pairs of vertices and their projections. Our goal is to find a plane distribution that minimizes $`E`$, thus optimizing the plane representation of the hypercube with respect to the distance between pairs. We have implemented a Montecarlo method to approach stochastically the optimal solution—the configuration of minimum “energy” $`E`$. Starting from a random initial configuration on the plane, each point performs a walk with fixed step length $`r`$ and directions chosen at random with uniform probability in $`[0,2\pi )`$. Each step of this walk produces a change in the configuration and, hence, in the distances $`d_{ij}`$, which implies a variation $`\mathrm{\Delta }E`$ in the energy. The new configuration is accepted with probability

$$p=\{\begin{array}{cc}\mathrm{exp}(\mathrm{\Delta }E/T)\hfill & \text{if }\mathrm{\Delta }E>0\hfill \\ 1\hfill & \text{otherwise},\hfill \end{array}$$
(2)
and rejected with probability $`1p`$. The “temperature” $`T`$ parametrizes this probability and allows the usual implementation of a simulated annealing, where the procedure starts with a high temperature that enables the system to explore a wide range in configuration space. Progressively, the temperature is reduced and the system freezes in one of the many local minima of the energy, typically not far away from the global minimum if the annealing is made slowly enough.

We have carried out the described procedure both interactively, reducing by hand the temperature while monitoring the configuration of the system on the computer screen, and automatically, by implementing a programmed reduction of the temperature. Our experiments show that essentially the same state is achieved in almost all the realizations. This implies that the energy landscape, while rugged, does not posses deep local minima that could capture the configuration far from the optimal one. The typical final configuration for $`L=10`$ ($`N=1024`$ points) is shown in Fig. 1(a). The self-similarity of its structure is remarkable, since no such property is present in the hypercube. Despite the appeal that this self-similar projection may have, it turns out that such projection is not well suited for our purpose. Vertices that are relatively near in the hypercube result rather far away in the projection. As an illustration, the first neighbors $`j`$ ($`h_{ij}=1`$) of a given vertex $`i`$ are shown in the figure. It is apparent that the Euclidean distance of some of them from the reference vertex is comparable with the size of the system. Moreover, many other vertices which should be farther away from vertex $`i`$, are in fact much closer.

One way to solve this difficulty is to modify the definition of $`E`$, such that near neighbors have more weight than distant neighbors. In fact, the energy (1) overemphasizes the effect of large distances. We have implemented the following simple alternative:

$$E=\underset{ij}{}\left(\frac{d_{ij}h_{ij}}{h_{ij}}\right)^2.$$
(3)
The final configuration, which we will term “homogeneous,” is shown in Fig. 1(b). Neighbor vertices of a vertex $`i`$ now result mapped onto points that surround the point $`i`$, which makes this projection much more satisfying. Certainly, however, some vertices result mapped near the border of the circle, and the arrangement of their neighbors is slightly different than that of vertices mapped in the middle of the set. We analyze below how this affects the distribution of distances.

A good characterization of the projections, and a quantitative way for comparing them, is the distribution of distances in each set. In the hypercube the distribution of distances to a vertex is the same for every vertex, and in fact is analytically found to be a binomial distribution. In the two-dimensional projections there is a different distribution for each point of the set. In Fig. 2(a) and (b) we show normalized distributions of distances for the self-similar and the homogeneous projections, respectively. In both figures, the black circles show the distribution of the distance to (any) vertex in the hypercube. Even though the distances form a discrete set, we show lines connecting the points to ease the reading of the graph. The other three curves shown in each plot correspond to the plane projections. The black squares correspond to an average on all the points in the sets. Triangles show averages performed on either the 10% of the points that form the external corona of the projection, or the 10% of its more central points. For the points of these subsets, still, all the distances to other points of the whole set are taken into account in the distributions.

The most immediate observation regarding Fig. 2 is the difference between the distributions in the two projections. The self-similar projection displays rugged distributions that reflect the hierarchical geometrical arrangement of the points. In the homogeneous projection, instead, the distributions are smooth, as in the hypercube. To this extent, the homogeneous projection can be said to represent more accurately the distribution of distances present in the hypercube. The distribution averaged over the whole set appears, however, skewed towards smaller distances, with a maximum around $`d=3`$, instead of the most represented distance $`h=5`$ of the hypercube. Interestingly, Fig. 2(b) shows that the *outer* 10% points considerably correct this skew. In other words, a point near the border of the circular array of the projection has a distribution of distances to the other points in the set which is rather similar to the distribution of a vertex of the hypercube.

An appraisal of the plane projections of the hypercube in a dynamical context results from the consideration of a diffusion process. Let us suppose that, at each time step, a random walker jumps from a vertex of the hypercube to one its neighbors with equal probability. The average distance $`D`$ from the initial site, as a function of time, is shown in Fig. 3(a) and (b) as black circles. The inset in both figures displays the same curve in double logarithmic scales, showing an initial behavior of the form $`D(t)t^{1/2}`$, like in a regular random walk in Euclidean space, followed by a saturation as the hypecube space is fully explored. The average distance as measured on the plane is shown in Figs. 3(a) and (b) for the self-similar and the homogeneous projections respectively. As expected from the distance distribution discussed above, the results for the plane projections depend on whether the initial point of the walker is at the border or at the center of the set. These two cases are shown in Figs. 3(a) and (b) as triangles pointing upward and downward, respectively. From this dynamical point of view, interior points behave equally bad in both projections. The most faithful representation of the process in a plane projection is the one given by one of the border points of the homogeneous set (Fig. 3(b), up triangles). Diffusion starting at these points behaves similarly as from points of the hypercube, both in the short and in the long time regimes, as seen in the linear and the logarithmic plots.

Our main goal of obtaining a sensible plane projection of the hypercube with the purpose of visualizing a dynamical process has been achieved, to an acceptable extent, by the homogeneous projection. Suppose that a dynamical phenomenon is taking place in a neighborhood of vertex $`P`$ of a hypercubical phase space. We need to build a homogeneous projection that maps vertex $`P`$ to a point at the border of the plane set. This is easily done by generating a projection at random and identifying one of the points at the border first. Suppose that one such point is $`Q`$. Then, each vertex $`I`$ of the hypercube is mapped to a point in the plane projection as

$$I(IP)Q,$$
(4)
where $``$ stands for the bitwise exclusive-or (XOR) operator. The projection obtained in this way provides a nice plane visualization substrate for the process.


text:
# Elemental abundances of metal poor carbon rich lead star: CS29497-030 Talk given at the meeting Early Galactic Chemical Evolution with UVES held at ESO Garching on 29-30th November 2002

T. Sivarani1 , P. Bonifacio1, P. Molaro1, R. Cayrel2, M. Spite 2, F. Spite2, B.Plez3, J. Andersen4, B. Barbuy5, T.C. Beers6, E. Depagne2, V.Hill2, P. François2, B. Nordström7,4, and F. Primas8
1 Istituto Nazionale di Astrofisicia - Osservatorio Astronomico di Trieste, Via Tiepolo 11, I-34131 Trieste, Italy
2Observatoire de Paris-Meudon, GEPI, F-92195 Meudon Cedex, France
3 GRAAL, Université de Montpellier II, F-34095 Montpellier Cedex 05, France
4 Astronomical Observatory, NBIfAFG, Juliane Maries Vej 30, DK-2100 Copenhagen, Denmark
5 IAG, Universidade de São Paulo, Departamento de Astronomia, CP 3386, 01060-970 São Paulo, Brazil
6 Department of Physics & Astronomy, Michigan State University, East Lansing, MI 48824, USA
7 Lund Observatory, Box 43, S-221 00 Lund, Sweden
8 European Southern Observatory (ESO), Karl-Schwarschild-Str. 2, D-85749 Garching b. München, Germany

We present here the abundance analysis of a metal poor carbon rich lead star, CS29497-030. High resolution and high signal to noise spectra were obtained using the UVES spectrograph on the 8.2m VLT-Kueyen telescope. The observations were made as a part of the Large Programme 165.N-0276, P.I. R. Cayrel. Abundance analysis was done using the latest version of the MARCS model atmospheres (Plez et. al. 1992) and the turbospectrum spectrum synthesis code. We have derived $`T_{\mathrm{eff}}`$ = 6650K from the FeI lines. Visible and infrared broad band colours using the Alonso et al. (1996) calibration, gives similar temperatures. A log $`g`$ value of 3.5 was obtained from the ionisation equilibrium of FeI and FeII, we remark that this gravity also satisfies the MgI/MgII, Ti I/Ti II and MnI/MnII equilibria, within errors. The abundance analysis indicates a metallicity, \[Fe/H\] = –2.7. A large overabundance of carbon (\[C/Fe\]=2.7) was found. We have also found large enhancement in the s-process elements and in particular lead shows an extremely high abundance of \[Pb/Fe\]=3.5, which makes this the star with the highest Pb/Fe ratio, up to date. The Pb/Ba ratio is found to be high (\[Pb/Ba\]=1.2) and the same is true for other second-peak s-process elements(e.g La, Ce, Nd). The star is a known spectroscopic binary with a period of 346 days (Preston & Sneden 2000). The abundance pattern suggests that CS 29497-30 has accreted matter from its companion, when it was in the AGB phase.


text:
# Casimir scaling hypothesis on the nonperturbative force in QCD vs. dual superconducting scenario of confinement11footnote 1 Talk given by Y. Koma at “PaNic02”, Osaka, Japan, Sep.30 - Oct.4, 2002

## 1 Introduction

From the SU(2) and SU(3) lattice QCD studies of the static potential between color charges in various dimensional representations, there arises the Casimir scaling hypothesis for the intermediate distance force . This hypothesis tells that the ratio of forces associated with various dimensional representations of color charges, for a given color group, is determined by that of eigenvalues of the quadratic Casimir operator for each representation.

It is expected that such a hypothesis is realized for the short distance force as described by one-gluon exchange, since the coupling is proportional to the quadratic Casimir operator. However, it is hard to imagine that this property is kept until intermediate distance where the nonperturbative effects set in. If the behavior of the ratio is governed exclusively by the group theoretical factor it should be manifest in arbitrary SU($`N`$) gauge theory. Recent studies of $`k`$-strings in SU(4) and SU(6) lattice gauge theories, however, do not support Casimir scaling . Therefore, it seems natural to consider that there are nonperturbative dynamics, rather than Casimir factor, which has inspired the Casimir scaling hypothesis for the intermediate distance force in early lattice investigations.

In this paper, we show that the dual superconducting picture of the nonperturbative QCD vacuum, as practically described by the dual Ginzburg-Landau (DGL) theory, provides us an understandable idea to explain the mechanism hidden behind the lattice data. In the dual superconducting vacuum, the color-electric flux is squeezed almost one dimensional tube, called the flux tube, due to the dual Meissner effect. For the quark and the antiquark system, then the flux tube is formed between the sources, leading to the linear potential, where its slope is identified as the string tension.

## 2 Flux-tube solution in the DGL theory

We examine the SU(3) case and calculate the string tensions of flux tubes in the DGL theory associated with static charges in various SU(3) representations systematically , starting from the manifestly Weyl symmetric form

$`_{\mathrm{DGL}}={\displaystyle \underset{i=1}{\overset{3}{}}}\left[{\displaystyle \frac{1}{6g^2}}{}_{}{}^{}F_{i\mu \nu }^{2}+\left|\left(_\mu +iB_{i\mu }\right)\chi _i\right|^2\lambda \left(\left|\chi _i\right|^2v^2\right)^2\right],`$ (1)
$`{}_{}{}^{}F_{i\mu \nu }^{}=_\mu B_{i\nu }_\nu B_{i\mu }+2\pi {\displaystyle \underset{j=1}{\overset{3}{}}}m_{ij}\mathrm{\Sigma }_{j\mu \nu }^{(e)}.`$ (2)
Here $`B_{i\mu }`$ and $`\chi _i`$ $`(i=1,2,3)`$ denote the dual gauge field and the complex scalar monopole field, respectively. The dual gauge fields within the Weyl symmetric expression have a constraint $`_{i=1}^3B_{i\mu }=0`$. The $`\mathrm{\Sigma }_{j\mu \nu }^{(e)}`$ $`(j=1,2,3)`$ in the dual field strength tensor is the color-electric Dirac string and their boundaries define the quark current $`j_{j\mu }^{(e)}`$ $`(j=1,2,3)`$.

The color-electric charge of the quark is specified by the weight vector of the SU(3) algebra, $`\stackrel{}{w}_j`$ ($`j=1,2,3`$), as $`\stackrel{}{Q}_j^{(e)}e\stackrel{}{w}_j`$, where $`\stackrel{}{w}_1=(1/2,\sqrt{3}/6)`$, $`\stackrel{}{w}_2=(1/2,\sqrt{3}/6)`$, and $`\stackrel{}{w}_3=(0,1/\sqrt{3})`$. On the other hand, the color-magnetic charges of the monopole fields $`\chi _i`$ are expressed by the root vectors of the SU(3) algebra, $`\stackrel{}{ϵ}_i`$, as $`\stackrel{}{Q}_i^{(m)}g\stackrel{}{ϵ}_i`$ ($`i=1,2,3`$), where $`\stackrel{}{ϵ}_1=(1/2,\sqrt{3}/2)`$, $`\stackrel{}{ϵ}_2=(1/2,\sqrt{3}/2)`$, and $`\stackrel{}{ϵ}_3=(1,0)`$. These color-electric and color-magnetic charges satisfy the extended Dirac quantization condition $`\stackrel{}{Q}_i^{(m)}\stackrel{}{Q}_j^{(e)}=2\pi m_{ij}`$, where $`m_{ij}=2\stackrel{}{ϵ}_i\stackrel{}{w}_j`$ and $`eg=4\pi `$. There are two mass scales in the DGL theory: the masses of the dual gauge boson $`m_B=\sqrt{3}gv`$ and the the monopole field $`m_\chi =2\sqrt{\lambda }v`$. In analogy to usual superconductors, their ratio, $`\kappa m_\chi /m_B`$, corresponds to the Ginzburg-Landau (GL) parameter, which describes the type of dual superconductivity of the vacuum.

The flux-tube solution is obtained by considering the cylindrical symmetric system with translational invariance along the $`z`$ axis. The fields are described as functions of radius $`r`$ and azimuthal angle $`\phi `$. Thus, we write the modulus of the monopole field as $`\varphi _i(r)=|\chi _i(r)|`$. The dual gauge field is parametrized as $`𝑩_i=[B_i^{\mathrm{reg}}(r)+B_i^{\mathrm{sing}}(r)]𝒆_\phi `$ where $`B_i^{\mathrm{reg}}(r)=[\stackrel{~}{B}_i^{\mathrm{reg}}(r)/r]`$ and $`B_i^{\mathrm{sing}}(r)=n_i^{(m)}/r`$ with $`n_i^{(m)}_{j=1}^3m_{ij}n_j^{(e)}`$. Here $`n_j^{(e)}`$ is the winding number of $`j`$-type color-electric Dirac string $`\mathrm{\Sigma }_{j\mu \nu }^{(e)}`$, which takes various integers depending on the representation of the SU(3) color group to which the charges belong. For a given representation $`D`$ of the SU(3) group, denoted by the Dynkin index $`[p,q]`$, we have $`\{n_1^{(e)},n_2^{(e)},n_3^{(e)}\}=\{p,q,0\}`$ .

The string tension of the flux tube is calculated as an energy per unit length in $`z`$ direction:

$`\sigma _D=2\pi {\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle _0^{\mathrm{}}}r𝑑r\left[{\displaystyle \frac{1}{3g^2}}\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{d\stackrel{~}{B}_i^{\mathrm{reg}}}{dr}}\right)^2+\left({\displaystyle \frac{d\varphi _i}{dr}}\right)^2+\left({\displaystyle \frac{\stackrel{~}{B}_i^{\mathrm{reg}}n_i^{(m)}}{r}}\right)^2\varphi _i^2+\lambda (\varphi _i^2v^2)^2\right].`$ (3)
In the Bogomol’nyi limit, $`\kappa =m_\chi /m_B=1`$, this expression is analytically evaluated as

$$\sigma _D=2\pi v^2\underset{i=1}{\overset{3}{}}\left|n_i^{(m)}\right|=4\pi v^2(p+q).$$
(4)
In this case, the ratio of the string tension between a higher and the fundamental representation $`[1,0]`$ is found to be $`d_D=\sigma _D/\sigma _\mathrm{𝟑}=p+q`$. In the general dual superconducting vacuum of type I ($`\kappa <1`$) and of type II ($`\kappa >1`$), one has to evaluate the whole expression (3) in its variational minimum by solving the field equations numerically.

In Fig. 1, we show the ratios of the string tensions of the flux tubes, $`d_D=\sigma _D/\sigma _\mathrm{𝟑}`$ for three values of the GL parameter, $`\kappa =1`$, $`3`$, and $`9`$ (numerically obtained for $`\kappa 1`$). We also plot the ratios of the string tensions obtained by the lattice simulations of Ref.  and the ratios of eigenvalues of the quadratic Casimir operator,

$`C^{(2)}(D)={\displaystyle \frac{1}{3}}(p^2+pq+q^2)+(p+q).`$ (5)
We find that the DGL result in the type II dual superconducting vacuum near $`\kappa =3`$ agrees well with all lattice data obtained in Ref. , albeit with big errors. The mechanism of the $`\kappa `$ dependence is understood as follows. In the Bogomol’nyi limit, $`\kappa =1`$, the ratio between the string tensions of a higher and the fundamental representation satisfies the flux counting rule: the string tension $`\sigma _D`$ is simply proportional to the number of the color-electric Dirac strings inside the flux tube, as seen from Eq. (4). With increasing $`\kappa `$, the interaction ranges of these fields get out of balance, and an excess of energy appears because of the interaction between fundamental flux tubes. This leads to systematic deviations from the counting rule. Note that the deviation of $`d_D`$ from the counting rule grows toward higher representations $`D`$, since the number of fundamental flux which coexist in the flux tube of representation $`D`$ increases as the sum $`p+q`$ of Dynkin indices. On the other hand, we also find that the DGL result at $`\kappa =9`$, for the deeply type II vacuum, uniformly reproduces Casimir-like ratios, through the deviations from the flux counting rule. In this analysis the higher dimensional flux tube in type II vacuum is assumed to be stable against splitting into fundamental ones. In principle, there must be a certain minimal $`q`$-$`\overline{q}`$ distance where such a effect is not negligible, depending on the GL parameter. However, in any case, the string tensions are saturated by the values at the Bogomol’nyi limit even if the splitting takes place.

## 3 Summary

We have studied the string tensions of flux tubes associated with static charges in various SU(3) representations in the DGL theory, based on a manifestly Weyl symmetric procedure. We have found that a GL-parameter near $`\kappa =3`$ reproduces the ratios of string tensions consistent with the lattice data . The DGL theory accidentally shows Casimir-like scaling for a deeply type II vacuum with $`\kappa =9`$. However, there is no direct relation to the eigenvalues of the Casimir operator. The mechanism of the systematic behavior of string tensions in the DGL theory can be understood as a result of the flux-tube dynamics. At present, it is not obvious that lattice data really contain such a dynamical effect. However, this example suggests that it is important to have more lattice results carefully interpreted without bias toward the Casimir scaling hypothesis.

## Acknowledgement

We are grateful to E.-M. Ilgenfritz and T. Suzuki, for useful discussions. This work is partially supported by the Ministry of Education, Science, Sports and Culture, Japan, Grant-in-Aid for Encouragement of Young Scientists (B), 14740161, 2002.


text:
# Cosmology with coalescing massive black holes

## Acknowledgments

We thank Shane Larson, David Merritt, Kristen Menou, and Steinn Sigurdsson for useful discussions. This work was support by NSF Grant PHY–9907949.

## References


text:
# Dynamics of a Bose-Einstein Condensate in an Anharmonic Trap

## References