text:
# 1 Introduction

## 1 Introduction

The drag coefficient for the turbulent boundary layer is determined as follows:

$$c_f^{}=\frac{2u_{}^2}{U^2},$$
$`[1]`$
where $`u_{}`$ is the friction velocity, $`u_{}=\sqrt{\tau /\rho }`$, $`\tau `$ the wall shear stress, $`\rho `$ the fluid density, $`U`$ is the free stream velocity. It is apparently the most practically important characteristic of the turbulent boundary layer. In the present note a correlation will be derived between the drag coefficient $`c_f^{}`$ and the effective Reynolds number $`Re`$, introduced previously (1,2) and determined as follows. It was shown (2) that when the free stream turbulence is low, the average velocity $`u`$ distribution in the intermediate region of the boundary layer between the viscous sublayer and the free stream consists of two different self-similar structures described by the scaling (power) laws $`\varphi =A\eta ^\alpha `$ and $`\varphi =B\eta ^\beta `$. Here $`\varphi =u/u_{}`$ and $`\eta =u_{}y/\nu `$ are dimensionless variables (recently the notations $`U^+`$ and $`y^+`$ instead of classical notations $`\varphi `$ and $`\eta `$ have become more popular), $`\nu `$ the fluid kinematic viscosity, $`y`$ the distance from the boundary. The interface between the two regions described by different scaling laws is sharp at $`y=\lambda `$, so that the viscous sublayer and the first region adjacent to it form the sharply bounded wall layer of thickness $`\lambda `$. Furthermore, we showed (2) that at sufficiently large Reynolds numbers the mean velocity distribution in the first region can be represented in the form obtained earlier for pipes (see (3))

$$\varphi =\frac{u}{u_{}}=\left(\frac{1}{\sqrt{3}}\mathrm{ln}Re+\frac{5}{2}\right)\left(\frac{u_{}y}{\nu }\right)^{\frac{3}{2\mathrm{ln}Re}}$$
$`[2]`$
where $`Re=U\mathrm{\Lambda }/\nu `$ is a Reynolds number, defined by a characteristic length scale $`\mathrm{\Lambda }`$ proportional to the thickness $`\lambda `$ of the wall layer (4). The procedure for determining the effective Reynolds number and the characteristic length $`\mathrm{\Lambda }`$, described in detail in (1,2), is not more complicated, perhaps even simpler than the procedure for determining the momentum thickness $`\theta `$ (1,4).

If the effective Reynolds number $`Re`$ is in fact a governing characteristic of the turbulent boundary layer, like the Reynolds number based on the average velocity and pipe diameter for flow in pipes, then the drag coefficient $`c_f^{}`$ should be a function of this parameter. Therefore it was natural to correlate the experimental data for drag coefficient and the effective Reynolds number determined for the same experiments. Such correlation will be presented and discussed below.

On the other hand, Watson, Hall and Anders (5) at the nasa Langley Research Center recently obtained a different correlation of the drag coefficient with the parameter traditional in turbulent boundary layer studies — the Reynolds number $`Re_\theta `$ based on the momentum thickness. A comparison of the two correlations is performed in the present Note.

## 2 Correlation of the drag coefficient and effective <br>Reynolds number for the turbulent boundary layer

According to the model presented by us previously (see (3)) the drag coefficient for the turbulent pipe flows $`c`$ as a function of the Reynolds number $`Re=\overline{u}d/\nu `$ is represented in the form

$$c=\frac{8u_{}^2}{\overline{u}^2}=\frac{8}{\mathrm{\Psi }^{2/(1+\alpha )}},\alpha =\frac{3}{2\mathrm{ln}Re}.$$

Here $`\overline{u}`$ is the mean velocity (bulk discharge per unit time divided by the cross-section area), and

$$\mathrm{\Psi }(\alpha )=\frac{e^{3/2}(\sqrt{3}+5\alpha )}{2^\alpha \alpha (1+\alpha )(2+\alpha )}.$$

As $`Re\mathrm{}`$, $`\mathrm{\Psi }(\alpha )(2e^{3/2}/\sqrt{3})\mathrm{ln}Re`$, so that asymptotically the relation is valid

$$c\frac{6}{e^3\mathrm{ln}^2Re}\frac{0.3}{\mathrm{ln}^2Re}.$$
$`[3]`$
The relation  is asymptotically covariant. This means that a replacement of $`Re`$ by certain $`Re^{}`$, equal to Constant times $`Re`$, leaves the asymptotics  invariant. Contrary to that, a redefinition of $`c`$, e.g. an introduction of a different factor before $`u_{}^2/\overline{u}^2`$ or the use of the maximum velocity instead of mean velocity $`\overline{u}`$, will change the Reynolds-number dependence of the drag coefficient.

Note that at large Reynolds number the expression for the drag coefficient should be identical for all shear flows provided the Reynolds number is properly determined. This identity was demonstrated for velocity profiles (1,2). Therefore it is logical to suggest that the expression for the drag coefficient in the boundary layer be of the form

$$c_f^{}=\frac{2u_{}^2}{U^2}=\frac{\text{Const}}{\mathrm{ln}^2Re}.$$
$`[4]`$
Here $`Re`$ is the effective Reynolds-number for the boundary layer, introduced in (1,2).

Processing of the experimental data by Naguib (6) and Nagib and Hites (see (7)) confirms (see Figure 1) the correlation , and gives the value of the Constant in  equal to 0.26, so that the correlation  takes the form

$$c_f^{}=\frac{0.26}{(\mathrm{ln}Re)^2}.$$
$`[5]`$
At the same time the data of Österlund (8) reveal a systematic deviation from (Figure 2). We suggest one of the reasons for this disagreement is the following fact: The friction velocity $`u_{}`$ was measured in (6,7) whereas it was calculated in (8) using the universal logarithmic law with inappropriately low values of the constants.

## 3 The correlation between the drag coefficient $`c_f^{}`$ and the Reynolds number $`Re_\theta `$ based on the momentum thickness

Watson, Hall and Anders (5) proposed recently the following correlation between the drag coefficient $`c_f^{}`$ and the traditional Reynolds number $`Re_\theta `$:

$$c_f^{}=0.0097Re_\theta ^{0.144}.$$
$`[6]`$
The correlation  was based on different experimental data covering the range of $`Re_\theta `$ from $`310^4`$ to $`610^5`$.

Our comparison showed (see Figures 3 and 4) that both the data of the experimentsof (6),(7) (Illinois Institute of Technology) and (8) (Roya l Institute of Technology, Stockholm) disagree systematically with the Langley correlation . The comparison of our correlation  with the results of Langley experiments was for us impossible due to the incompleteness of the data available to us.

We notice that the power type asymptotics with small values of the exponent similar to  can be approximated by formulas similar to  and vice versa. Indeed, let us approximate $`1/\mathrm{ln}^2x`$ by a power function

$$\frac{1}{\mathrm{ln}^2x}=Gx^\gamma .$$

The constants $`G`$ and $`\gamma `$ can be determined for instance from the condition that at a certain point $`x_0`$ the functions $`Gx^\gamma `$ and $`1/\mathrm{ln}^2x`$ be equal as well as their derivatives. This condition gives $`\gamma =2/\mathrm{ln}x_0`$, and $`G=x_0^\gamma (\mathrm{ln}x_0)^2`$. So, if the correlation  has a physical meaning, the empirically obtained correlations of the power type similar to  can be in fact their approximations.

## 4 Conclusion

A correlation between the drag coefficient and the effective Reynolds number for the turbulent boundary layers at large Reynolds numbers is proposed. It is in satisfactory agreement with the data of Naguib (6), and Nagib and Hites (7).

The data of Österlund (8) for the drag coefficient as a function of the effective Reynolds number are in systematic disagreement with the correlation . Therefore there is a systematic disagreement between the data (6),(7) on one hand, and (8) on the other.

All the data of (6), (7) and (8) are is substantial systematic disagreement with the correlation  proposed in (5). A detailed analysis of the experimental data of (5) is needed. If such a correlation can be firmly established, it will provide the relation between the characteristic length scale $`\mathrm{\Lambda }`$ and the momentum thickness $`\theta `$.

Acknowledgments. This work was supported in part by the Applied Mathematics subprogram of the U.S. Department of Energy under contract DE–AC03–76–SF00098, and in part by the National Science Foundation under grants DMS 94–16431 and DMS 97–32710.

References

1. Barenblatt, G.I., Chorin, A.J., Hald, O. and Prostokishin, V.M. (1997). Proc. Nat. Acad. Sci. USA 94, 7817–7819.
2. Barenblatt, G.I., Chorin, A.J. and Prostokishin, V.M. (2000). J. Fluid Mech. 410, 263–283.
3. Barenblatt, G.I., Chorin, A.J. and Prostokishin, V.M. (1997). Appl. Mech. Rev. 50, 413–429.
4. Barenblatt, G.I., Chorin, A.J. and Prostokishin, V.M. (2000). Proc. Nat. Acad. Sci. USA 97, 3799–3802.
5. Watson, R.D., Hall, R.M. and Anders, J.B. (2000). Review of Skin Friction Measurements Including Recent High-Reynolds Number Resutls from nasa Langley NTF. Paper AIAA–20000-2392, 1–20.
6. Naguib, A.M. (1992). Inner- and outer-layer effects on the dynamics of a turbulent boundary layer, Ph.D. Thesis, Illinois Inst. of Technology, Chicago, Ill.
7. Hites, M.H. (1997). Scaling of High-Reynolds Number Turbulent Boundary Layers in the National Diagnostic Facility. Ph.D. Thesis, Illinois Inst. of Technology, Chicago, Ill.
8. Österlund, J.M. (1999). Experimental studies of zero pressure-gradient turbulent boundary layer flow, Ph.D. Thesis, Royal Inst. of Technology, Stockholm, Sweden.

> Figure Captions
> Figure 1: The correlation  is in general agreement with the experimental data of Naguib (6) (diamonds) and Nagib and Hites (7) (squares).
>
> Figure 2: There is a systematic substantial disagreement of the correlation  and the data of Österlund (8) (circles).
>
> Figure 3: There is a systematic disagreement between the correlation  and the experimental data (6,7).
>
> Figure 4: There is a systematic disagreement between the correlation  and the experimental data (8).


format: text

text:
# Hartree-Fock-Bogoliubov Calculations of the Rotational Band of the Very Heavy 254No Nucleus

11institutetext: Centre d’Études Nucléaires de Bordeaux Gradignan, CNRS-IN2P3 and Université Bordeaux-I, BP 120, F-33175 Gradignan, France 22institutetext: Institute of Nuclear Research and Nuclear Energy (Bulgarian Academy of Sciences), Tzarigradsko Chaussee 72, 1784 Sofia, Bulgaria

## Abstract

We report on Hartree-Fock-Bogoliubov (HFB) calculations of the ground-state rotational band of the heavy nucleus 254No recently observed experimentally. The calculated quadrupole deformation is consistent with the experimental value of $`\beta =0.27`$ and is almost constant over the whole band. We also reproduce fairly well the excitation spectra and moments of inertia of this isotope up to the maximal experimentally observed state of spin 20. The rather high stability of this nucleus against fission is illustrated by the deformation energy curve providing very high fission barriers at zero spin within the HFB and HFB plus Lipkin-Nogami formalisms. The variation of these barriers with increased angular velocities is also studied.

dedicatory: This paper is dedicated to the memory of Jordanka Piperova

After the Cohen, Plasil and Swiatecki CPS  seminal study of the stability of rotating liquid drops, very heavy nuclei with mass greater than 250 are generally expected to survive rotational excitation only up to rather low spins. Of course, large shell-correction energies added to the bulk liquid drop estimates may lead to sizable enhancements of the stability against fission in some cases. Observation of relatively high-spin states in very heavy nuclei thus provides a very useful tool to assess any theoretical description of shell effects in this mass region, which is of primary importance for a good prediction of shell-stabilized superheavy elements.

The experimental observation in the ground-state band of the 254No nucleus of rotational states from spins 4 to 18 $`\mathrm{}`$ (and tentatively 20 $`\mathrm{}`$) has been reported recently in a series of papers Expe:I ; Expe:II ; Expe:III . These experiments have been performed using 48Ca beams from 130 to 219 MeV on 208Pb targets. Two of them (including the more recent) were conducted at the Argonne National Laboratory with the Gammasphere $`4\pi `$ Ge detector array, while the other made use of four Ge clover detectors and was performed at the University of Jyväskylä. In all three experiments, the 254No nuclei were implanted into a position-sensitive Si strip-detector placed at the focal plane. The $`\gamma `$-rays coming from the 254No were unambiguously identified from the coincidence with the $`\alpha `$-decay chain in the Si detector. As the transition to the ground-state have not been observed, the spin assignment relies on the use of a Harris parameterization for both $`J^{(1)}`$ and $`J^{(2)}`$ moments of inertia.

In ref. Expe:III , the authors proposed a new method to deduce the fission barrier height $`B_f`$ based on a reconstruction of the entry distribution of the evaporation residues in spin and excitation energy. Since the fission process above the saddle-point energy should be much favored compared to the $`\gamma `$-emission, they claim that for any given observed spin, the maximal excitation energy of the entry distribution should lie below the saddle-point energy. For each spin the end-point of the excitation energy distribution gives thus a lower bound for the fission barrier height which is estimated to be greater than 5 MeV around $`I=12`$ $`\mathrm{}`$.

In a recent paper involving some of the authors Ldm:rot , it has been shown using state-of-the-art semiclassical calculations, namely in the extended Thomas-Fermi (ETF) approach, that the rotating liquid drop model had to be refined by including a spin-dependence of its parameters which, as a result, provides a fission stability of heavy nuclei upon increasing the spin which is enhanced with respect to what is stated in ref. CPS . It is the aim of the present paper to go beyond this work, describing altogether bulk properties as well as shell and pairing effects within the fully self-consistent microscopic Hartree-Fock-Bogoliubov (HFB) formalism, in the particular case of 254No.

In our calculations, we have used the triaxial Hartree-Fock-Bogoliubov code developed by Laftchiev et al. HFB:code ; Hristo:these . This code is an extension of the Hartree-Fock Routhian code presented in refs. HF:code ; Period  which uses Skyrme-type effective interactions and assumes parity and signature as symmetry operators for the one-body Hamiltonian. There, solutions breaking time-reversal symmetry are described through a decomposition of the single-particle wavefunctions on an axially symmetric harmonic oscillator basis. The triaxial character inherent to the solutions of this formalism is taken care of by a decomposition of the various densities as Fourier series in the azimuthal angle. This particular choice provides shorter computation-times than usual triaxial codes as for instance those of refs. Dudu:I ; Dudu:II . The pairing correlations are described à la HFB with particle-number projection tentatively taken into account within the approximate Lipkin-Nogami (LN) scheme Nogami ; Kamlah . The latter method which is widely used for instance in the context of fully self-consistent microscopic calculations (see e.g., ref. BCSLN  for time-reversal symmetric Hamiltonians, or ref. Gall  for Routhian calculations) results in the addition of a constraint on the second-order fluctuations of the particle-number. However it has been shown Peru  in the $`A190`$ mass region from calculations using the Gogny force that this approximation could produce somewhat inconsistent results.

The HFB equations may be cast into a usual eigenvalue form in a doubled space where the $`2\times 2`$ “Hamiltonian” is defined in terms of the Hartree-Fock Hamiltonian $`h`$ and of the pairing potential $`\mathrm{\Delta }`$, whereas the eigensolutions are the usual $`U`$ and $`V`$ matrices (see for the notations, e.g, RingSchuck ). Similarly to what was done in ref. Gall , these equations are solved in two steps at each iteration. First we solve simple Hartree-Fock equations (that is the HFB equations with vanishing pairing field) to obtain the eigenstates of the Hartree-Fock–like Hamiltonian. Then the original equations are solved in a truncated configuration space (through an energy cutoff) of the Hamiltonian eigenstates previously determined.

For the particle-hole channel, we have used the SkM parameterization of the Skyrme interaction Skmstar  since it was originally fitted to provide a good description of the 240Pu fission barrier and is thus rather well adapted to describe the stability against fission in very heavy nuclei. In the particle-particle channel, we have simply used the “modified seniority pairing force” previously described by Gall et al. in ref. Gall . At the present stage of our investigations of rotational properties in very heavy nuclei, we have not attempted to make a global fit over a variety of different nuclei to provide a pairing force which could be deemed to be somewhat universal in this mass region. We rather satisfied ourselves with a parameterization yielding the right moment of inertia at very low spins to assess the validity of our approach in reproducing the transition energies (and moment of inertia at higher spins). The matrix elements $`G_n`$ and $`G_p`$ for neutrons and protons are defined with the usual prescription $`G_q=g_q/(11+N_q)`$, $`N_q`$ being the number of particles in the charge state $`q`$. The value of $`g_n`$ ($`g_p`$ resp.) is 14.3 MeV (15.5 MeV resp.) in the pure HFB case and 10.65 MeV (14.1 MeV resp.) in the HFB+LN case with an energy cutoff 6.2 MeV above the Fermi level in both cases. No smeared boundary conditions for this configuration space have been taken into account yet which may locally create some minor technical deficiencies.

The calculations have been performed using the usual deformation-dependent truncation scheme equivalent to a 13 major shells spherical basis, corresponding to more than 450 harmonic oscillator states (for both signatures). The basis deformation parameters that have been used in our HFB calculations were obtained through a minimization of the energy calculated within an axially and time-reversal symmetric HF+BCS formalism. At zero angular velocity non axially-symmetrical solutions have been forcefully searched out (due to the symmetry properties of the interaction, one would never explore as well known non-axial solutions starting from a first iteration axially-symmetrical ansatz). At finite spins, the solutions explicitly break the axial symmetry. As a result however, the solutions are almost axially symmetric since the $`\gamma `$ angle never exceeds 1 along the rotational band which allows us to only take into account the 5 lowest-order Fourier series in the densities decomposition (see ref. HF:code ).

Recently, T. Duguet et al. Nobel:bonche  have presented the results of HFB+LN calculations which could be considered to be in essence rather similar to ours with two practical differences which are however not very substantial provided that the technical work is performed adequately \[(i) solutions of the HFB equations are computed on a spatial grid, (ii) a zero-range force is used in the particle-particle channel\]. However they used for the particle-hole channel the SLy4 parameterization Sly4a ; Sly4b  which is known to give in average (i.e., semiclassically) significantly too high fission barriers PrivCom . This may be considered as a drawback to study fission stability properties. Another series of calculations for the same nucleus has been published quite recently Egido  making use of the D1S Gogny force parameterization Gogny:force  which is known to be well suited to study fission barriers. No corrections for the particle number symmetry violation, e.g., of the Lipkin-Nogami type, have been taken into account.

Along the rotational band of 254No, the mass quadrupole moment $`Q_0`$ varies very slightly between 32.9 and 32.7 b in our calculations. This corresponds to a deformation parameter of $`\beta =0.27`$ identical to the experimental value. It is to be noted that we are there in perfect agreement with the experimental results Expe:I ; Expe:II ; Expe:III  as well as with the theoretical results of T. Duguet et al. Nobel:bonche  yielding $`Q_{20}=32.8`$ barn. The calculations of Egido and Robledo on the other hand yield slightly higher value of $`Q_{20}`$, resulting in a $`\beta `$ value close to 0.29.

We present in fig. 1 the evolution of the pairing energies, defined as a quantity proportional to the trace of the product of the abnormal density $`\kappa `$ with the pairing potential $`\mathrm{\Delta }`$, for protons and neutrons with respect to the angular velocity $`\mathrm{\Omega }`$. It is well known from Nilsson-type calculations that for a deformation parameter of $`\beta 0.3`$ the $`N=152`$ neutron number is magic Magic . This prediction is substantiated in our case by the rather low value of the neutron pairing energy compared to the proton corresponding energy. As usual in HFB calculations, the pairing energy is seen to decrease with increasing angular velocity. In particular, the neutron pairing energy vanishes above $`\mathrm{}\mathrm{\Omega }=0.22`$ MeV. Within the same angular velocity range, the proton pairing energy is reduced by more than 40 percent.

It is clear from fig. 1, that around and above rotational frequencies of the order of $`\mathrm{}\mathrm{\Omega }=0.1`$ MeV the results of pure HFB calculations are rather dubious as far as the neutron pairing correlations are concerned. This remark is clearly applicable also to the calculations of Egido and Robledo. It is generally expected that the Lipkin-Nogami corrected HFB calculations should be better adapted to such a situation of weak pairing correlations. However, as already mentioned, the calculations of Peru et al. (See, e.g., Peru ) in the superdeformed $`A190`$ region show that the quality of the results of such an approach is rather unpredictable. The application to such a Routhian approach of pairing calculations conserving explicitly the particle number, e.g., as those of Pillet et al. Pillet ; Pillet:phd  would be of great interest there. In view of this expected inadequacy of pure HFB approach we have therefore limited the discussion below of calculated moments of inertia to HFB+LN results.

The kinematic ($`J^{(1)}`$) and dynamic ($`J^{(2)}`$) moments of inertia are obtained from our calculations using the formulae

$`J^{(1)}`$ $`={\displaystyle \frac{\widehat{j}_x}{\mathrm{\Omega }}}`$ (1)
$`J^{(2)}`$ $`={\displaystyle \frac{\widehat{j}_x}{\mathrm{\Omega }}}`$ (2)
where $`\widehat{j}_x`$ is the total angular momentum operator and $`\mathrm{\Omega }`$ is the angular velocity. The derivative involved in the $`J^{(2)}`$ calculations have been taken over two consecutive even-spins. Both are plotted in fig. 2 as functions of the angular velocity together with their experimental counterparts.

The kinematic moment of inertia $`J^{(1)}`$ is, in fact, the expression of the linear response of the system under the constraint of some collective variable. Nevertheless, we have actually assumed that it could be computed from the wavefunctions yielding the corrected energy (i.e., bypassing an explicit treatment of LN corrections specific to $`j_x`$ in the spirit of ref. BCSLN ). It may also be argued that due to the dependence of $`\lambda _2`$ upon the angular velocity $`\omega `$, equation (2) is not correct and dynamic moments of inertia should be evaluated through second derivatives of the energy with respect to the spin. However, we have checked that this $`\omega `$-dependence is very small and can be safely neglected (at least up to about $`\omega =0.1`$ MeV).

It is seen on fig. 2 that we are reproducing rather well the experimental trend for both the $`J^{(1)}`$ and $`J^{(2)}`$ moments of inertia. The agreement of our results is better than the one obtained by Egido and Robledo for the only moment which they report, namely $`J^{(1)}`$. It is also better than the only moment reported by Duguet et al., namely $`J^{(2)}`$ in this case. The latter is quite remarkable because, as it can be seen on fig. 3, our single particle spectra are somewhat different for the neutron states. Whereas these authors find important deformed gaps feature around $`N=150`$ and $`N=152`$ (see fig. 5 in Nobel:bonche ), we get only a comparable $`N=152`$ gap plus a $`N=170`$ gap at low spins which is not present in the paper of Duguet et al..

On fig. 4 we display the fission barriers obtained at zero total angular momentum within the HFB and HFB+LN approaches allowing triaxial deformation. The former yields a first fission barrier which is slightly lower than what is obtained for the latter. Comparing HFB+LN results now for two different total angular momenta $`I=0`$ $`\mathrm{}`$ and $`I=12`$ $`\mathrm{}`$ we find that the first fission barrier height is slightly increased by 500 keV where one would expect it to be decreased. This phenomenon is due to a higher single-particle level density at the top of the first barrier for $`I=12`$ $`\mathrm{}`$ as compared to what is obtained at zero-spin. At such low spins, shell structure effects are thus able to mask the anti-binding effect of the rotation leading to an extra fission instability With the basis size in use in our calculations we are not a priori able to provide reliable relative energies in the second fission barrier region, as noted years ago Flocard:basis .

To conclude it appears that the SkM parameterization of the Skyrme interaction (complemented by a seniority force in the p-p h-h channel) is able to provide a good reproduction of the spectroscopic data in the ground state rotational band of the 254No nucleus within the HFB plus Lipkin-Nogami approach. This result is all the more remarkable that similar calculations Nobel:bonche  using a different effective force and yielding a rather different shell structure in the relevant deformation range, lead to a similarly (may be slightly less) good reproduction of the experimental data. Even though it is true that reproducing such collective rotation data provides a much wanted test of effective forces to be used when predicting the stability properties of neighbouring superheavy nuclei, such a convergence of results from seemingly rather different effective forces should mitigate the hope to get so far here a completely stringent benchmark.

The calculations of ref. Egido  have studied in detail the fission properties of this nucleus, namely they underline its remarkable, and somewhat unexpected, fission stability which could partly be explained by a change of saturation properties of nuclear matter in presence of a centrifugal field as advocated in ref. Ldm:rot . Our results, limited in both spins and deformations as compared to those of ref. Egido  yield similar conclusions on this stability. A pending point common to all these calculations (those of refs. Nobel:bonche ; Egido  as well as ours) is the pairing correlation treatment in weak pairing regions as encountered upon increasing the angular velocity. Surely, treatment of such correlations within an approach that explicitly conserves the particle number, as the one proposed in refs.Pillet ; Pillet:phd , would greatly improve the validity of the results.

However whatever their mid stream character, the rather satisfactory reproduction of the difficult, and for that reason somewhat scarce, experimental results so far obtained by self consistent calculations as those of refs. Nobel:bonche ; Egido  and the one presented here, should constitute an incentive to experimentalists to provide more data so as to improve our knowledge of the effective interaction to be used. This certainly constitutes a very timely endeavour.


text:
# Characterization of the probabilistic traveling salesman problem

## I Introduction

Many real systems present problems of stochastic optimization. These include communications networks, protein design fink1  and oil field models jonsbraten , in all of which uncertainty plays a central role. We will consider the case where the outcome $`g(x,\omega )`$ depends not only on parameters $`x`$ to be chosen, but also on unknowns $`\omega `$. We can only average with respect to these unknowns, aiming to find the ‘solution’ $`x`$ which optimises the average outcome. Thus we seek to find $`xX`$ which minimizes

$$\overline{g}(x)=g(x,\omega )f(\omega )𝑑\omega $$
(1)
where $`X`$ is the solution space of the problem and $`f(\omega )`$ is the probability distribution of the uncertain variables.

Stochastic optimization was borne out of an idea by Robbins & Monro robbins . They considered solving the problem of finding

$$G(x)=\alpha $$
(2)
where $`G`$ is some monotonic function of $`x`$ and $`\alpha `$ is a parameter. $`G`$ is not known directly, but can only be estimated. Their technique to solve this problem is called stochastic approximation, and a number of variants of this scheme have since been developed benveniste ; kushner1 ; l'ecuyer ; kushner2 .

Stochastic optimization is a more general situation since the function to be minimized may have many local minima. We may classify the techniques to solve stochastic optimization problems into two classes - exact methods and heuristics. Heuristics are more appropriate to NP-complete problems, and these are the problems on which we focus in this paper. A number of heuristics already exist to tackle stochastic optimization problems gong ; yan ; devroye ; yakowitz ; andradottir . Many of these are developments from simulated annealing haddock ; bulgak ; alrefaei ; alkhamis , which has been shown gutjahr  to solve stochastic optimization problems with probability 1, provided $`\overline{g}(x)`$ can be estimated with precision greater than $`O(t^\gamma )`$ for time step $`t`$, where $`\gamma >1`$. A number of authors bulgak ; alrefaei ; alkhamis ; gelfand  have used a modified simulated annealing algorithm in which the acceptance probability is modified to take some account of the precision of the estimates of $`\overline{g}(x)`$, and in these cases there are a number of convergence results gelfand ; alkhamis .

Stochastic annealing fink1  is a modified simulated annealing algorithm which differs from the above approaches in two key ways. Firstly the noise present in estimates is positively exploited as mimicking thermal noise in a slow cooling, as opposed to being regarded as something whose influence should be minimised from the outset. Secondly, stochastic annealing can be modified to give exact simulation of a thermal system. Although this is not specifically ruled out by the earlier approaches, no attempt has been made with them to satisfy this condition.

In stochastic annealing we estimate $`\overline{g}(x)`$ by taking $`r`$ repeated, statistically independent, measurements of $`g(x,\omega )`$ each of which we call an instance. For the implementation of stochastic annealing used in this paper we accept all moves for which one estimate (based on $`r`$ instances) for a new state is more favourable than an equivalent estimate for the old. This simple procedure does not exactly simulate a thermal system, where the acceptance probabilities should obey

$$\frac{P_{AB}}{P_{BA}}=e^{\beta \mathrm{\Delta }\mu }$$
(3)
where $`\beta =\frac{1}{k_BT}`$ and $`\mathrm{\Delta }\mu `$ is the exact difference in $`\overline{g}(x)`$ between states $`A`$ and $`B`$. However if we assume that our estimate change is Gaussian distributed around $`\overline{g}(x)`$ with standard deviation $`\frac{\sigma }{\sqrt{r}}`$ , where $`r`$ is the number of instances used for each estimate, then it follows that the acceptance probability is fink1

$$P_{AB}^G=\frac{1}{2}\left[1\text{erf}\left(\frac{\sqrt{r}\mathrm{\Delta }\mu }{\sqrt{2}\sigma }\right)\right].$$
(4)
The approximation to a thermal acceptance rule is then quite good since

$$\begin{array}{ccc}ln\left(\frac{P_{AB}}{P_{BA}}\right)\hfill & =ln\left(\frac{1\text{erf}\left(\frac{\sqrt{r}\mathrm{\Delta }\mu }{\sqrt{2}\sigma }\right)}{1+\text{erf}\left(\frac{\sqrt{r}\mathrm{\Delta }\mu }{\sqrt{2}\sigma }\right)}\right)\hfill & \\ & \beta _G\mathrm{\Delta }\mu \frac{4\pi }{48}(\beta _G\mathrm{\Delta }\mu )^3\mathrm{}\hfill & \end{array}$$
(5)
where

$$\beta _G=\frac{\sqrt{8r}}{\sqrt{\pi }\sigma }$$
(6)
identifies the equivalent effective temperature. The small coefficient ($`0.02`$) of the cubic term in eq. 5 makes this a rather good approximation to true thermal selection.

Increasing sample size $`r`$ means that we are more stringent about not accepting moves that are unfavourable, equivalent to lowering the temperature, which is quantified by eq. 6 for the Gaussian case. As with standard simulated annealing rees ; nulton ; salamon1 , the question of precisely what cooling schedule to use remains something of an art.

## II Probabilistic Traveling Salesman Problem (PTSP)

We adopt the PTSP as a good test-bed amongst stochastic optimization problems, in much the same way as the TSP has been considered a standard amongst deterministic optimization problems. The PTSP falls into the class of NP-complete problems bertsimas4 , and the TSP is a subset of the PTSP.

The original traveling salesman problem (TSP) is to find the shortest tour around $`n`$ cities, in which each city is visited once. For small numbers of cities this is an easy task, but the problem is NP-complete: it is believed for large $`n`$ that there is no algorithm which can solve the problem in a time polynomial in $`n`$. Consideration of the traveling salesman problem began with Beardwood *et al.* beardwood . They showed that in the limit of large numbers of cities which are randomly distributed on the unit square, the optimal tour length ($`L_{\text{TSP}}`$) follows steele1

$$E(L_{\text{TSP}})=\beta _{\text{TSP}}\sqrt{n}+\alpha $$
(7)
where $`\beta _{\text{TSP}}`$ and $`\alpha `$ are constants. Here and below $`E(L)`$ denotes the quantity $`L`$ averaged, after optimization, with respect to different city positions, randomly placed on the unit square. Numerical simulation lee  gives $`\beta _{\text{TSP}}=0.7211(3)`$ and $`\alpha =0.604(5)`$ as estimates when $`n50`$. Significant divergence from this behaviour is found for $`n10`$, but numerical estimates can be found quickly (see appendix).

The probabilistic traveling salesman problem (PTSP), introduced by Jaillet jaillet1 ; jaillet2 , is an extension of the traveling salesman problem to optimization in the face of unknown data. Whereas all of the cities in the TSP must be visited once, in the PTSP each city only needs to be visited with some probability, $`p`$. One first decides upon the order in which the cities are to be visited, the ‘a priori’ tour. Subsequently, it is revealed which cities need to be visited, and those which do not need to be visited are skipped to leave a ‘pruned tour’. The order in which the cities are to be visited is preserved when pruning superfluous cities. The objective is to chose an a priori tour which minimizes the average length of the pruned tour. It is clear from figure 1 that near optimal a priori tours may appear very different for different values of $`p`$.

In our terminology, the average pruned tour length is averaged over all possible instances of which cities require to be visited. This was given by Jaillet as jaillet1

$$\overline{L}_{\text{pruned}}=\underset{q=0}{\overset{n2}{}}p^2(1p)^qL_t^{(q)}$$
(8)
where

$$L_t^{(q)}=\underset{j=1}{\overset{n}{}}d(j,1+(j+q)_{\text{mod}n})$$
(9)
is the sum of the distances between each city and its $`(q+1)^{th}`$ following city on the a priori tour, and the factors $`p^2(1p)^q`$ in the preceding equation simply give the probability that any particular span skipping $`q`$ cities occurs in the pruned tour. Jaillet’s closed form expression for the average pruned tour length renders the PTSP to some extent accessible as a standard (but still NP complete) optimization problem, and provides some check on the PTSP results by stochastic optimization methods.

It has been conjectured bertsimas4  that, in the limit of large $`n`$, the PTSP strategy is as good as constructing a TSP tour on the cities requiring a visit, the re-optimization strategy. This would mean that

$$\underset{n\mathrm{}}{lim}\left(\frac{E\left(\overline{L}_{\text{pruned}}\right)}{\sqrt{np}}\right)=\beta _{\text{TSP}},$$
(10)
where $`E(\overline{L}_{\text{pruned}})`$ is the pruned tour length further averaged over city positions after optimisation, which we will refer to as the expected pruned tour length. Figure 2 shows the expected pruned tour length divided by the expected re-optimized tour length. Since this quantity is tending towards a value significantly greater than 1 for $`p<1`$ it demonstrates that the PTSP strategy can be worse than the re-optimization strategy. Jaillet jaillet1  and Bertsimas *et al.* bertsimas2  have also shown that there is a limit to how much worse it can be, with

$$\underset{n\mathrm{}}{lim}\left(\frac{E\left(\overline{L}_{\text{pruned}}\right)}{\sqrt{np}}\right)=\beta _{\text{pruned}}(p)$$
(11)
where

$$\beta _{\text{TSP}}\beta _{\text{pruned}}(p)Min(0.9212,\frac{\beta _{\text{TSP}}}{\sqrt{p}}).$$
(12)
One attempt to solve the PTSP using an exact method was taken by Laporte *et al.* laporte  who introduced the use of integer linear stochastic programming. Although use of algorithms which may exactly solve the PTSP are useful, they are always very limited in the size of problem which may be attempted. Furthermore, the stochastic programming algorithm even fails to solve the PTSP on certain occasions, thus the accuracy of any statistics that would be generated using this method is dubious.

Three studies have used heuristics to solve the PTSP bertsimas1 ; rossi ; bertsimas5 . None of these studies used global search heuristics, and all were very restricted in the problem size attempted due to computational cost. The evaluation of a move for the PTSP using equation 8 involves the computation of $`O(n^2)`$ terms compared to $`O(1)`$ computations to evaluate a move in the TSP. Thus, to solve a 100 city problem for the PTSP would take $`O(10,000)`$ times longer than it would to solve a 100 city problem for the TSP. It should be noted, however, that it is only possible to make this comparison due to the relative simplicity of the PTSP. For many more stochastic optimization problems, standard optimization techniques are simply not applicable.

## III Form Of The Optimal Tour & Scaling Arguments

Optimal a priori PTSP tours for small $`p`$, as exemplified in figure 1 for $`p=0.1`$, resemble an “angular sort” - where cities are ordered by their angle with respect to the centre of the square. Bertsimas bertsimas1  proposed that an angular sort be optimal as $`p0`$, but we can show this to be false by comparison to a space-filling curve algorithm which is generally superior as $`n\mathrm{}`$. Such an algorithm was introduced by Bartholdi *et al.* bartholdi  using a technique based on a Sierpinski curve.

For the angular sort with $`np1,`$ the probability of two cities being nearest neighbours on the pruned tour will be vanishingly small for cities which are separated from each other by a large angle on the a priori tour. This means that only cities that are separated by a small angle contribute significantly to eq. 8. Thus for an $`n`$ city tour chosen by angular sort, we may approximate

$$L_t^{(q)}L_on$$
(13)
where $`L_o`$ is some fraction of the side of the unit square, since cities which are sorted with respect to angle will be unsorted with respect to radial distance. This leads to

$$E(\overline{L}_{\text{ang}})L_onp^2\underset{q=0}{\overset{n2}{}}(1p)^q.$$
(14)
For $`np1`$ and $`p1`$, we then find that the angular sort yields

$$E(\overline{L}_{\text{ang}})L_onp.$$
(15)
By contrast it has been shownbertsimas2  that

$$\frac{E(\overline{L}_{\tau _{\text{sf}}})}{E(\overline{L}_{\text{Reopt}})}=C$$
(16)
with probability 1, where $`E(\overline{L}_{\tau _{\text{sf}}})`$ is the expected length of a tour generated by a heuristic based on the Sierpinski curve and $`E(\overline{L}_{\text{Reopt}})`$ is the expected length for the re-optimization strategy. Using previous computational results bertsimas2 ; lee , we estimate $`C1.33`$, which is worse than we achieve using stochastic annealing. Hence, $`E(\overline{L}_{\tau _{\text{sf}}})`$ is given by

$$E(\overline{L}_{\tau _{\text{sf}}})=O(\sqrt{np})$$
(17)
which leads to

$$\frac{E(\overline{L}_{\tau _{\text{sf}}})}{E(\overline{L}_{\text{ang}})}=O(\frac{1}{\sqrt{np}}).$$
(18)
So for large enough $`np`$, the angular sort is not optimal.

From inspection of near-optimal PTSP tours such as fig. 1, we propose that the tour behaves differently on different length scales; the tour being TSP-like at larger length scales, but resembling a locally directed sort at smaller length scales. We may construct such a tour and use scaling arguments to analyse both the pruned and a priori lengths of the optimal tour. Consider dividing the unit square into a series of ‘blobs’, each blob containing $`1/p`$ cities so that of order one city requires a visit. The number of such blobs is given by

$$Nnp$$
(19)
and for these to approximately cover the unit square their typical linear dimension $`\xi `$ must obey

$$N\xi ^21.$$
(20)
Since each blob is visited of order once by a pruned tour, we can estimate the expected pruned tour length to be

$$E(\overline{L}_{\text{pruned}})N\xi \sqrt{np}$$
(21)
which we will see below is verified numerically. We can similarly estimate the a priori tour length to be $`n`$ times the distance between two cities in the same blob. Thus, the expected a priori tour length is

$$E\left(L_{\text{a priori}}\right)n\xi \sqrt{\frac{n}{p}}$$
(22)
which is more difficult to confirm numerically.

## IV Computational Results For The PTSP

We have investigated near optimal PTSP tours for a range of different numbers of cities, and various values of $`p`$. We used stochastic annealing with effective temperatures in the range $`kT=0.070.01`$, corresponding to sample sizes in the range $`r=2500`$. Between $`10`$ and $`80`$ different random city configurations were optimized ($`80`$ configurations of 30 cities, $`40`$ configurations of 60 cities, $`20`$ configurations of 90 cities and $`10`$ configurations for $`n120`$ cities).

Figure 3 shows a master curve for the expected pruned tour length divided by $`\sqrt{np}`$. The shift factors have a linear fit and the data are consistent with

$$\frac{E(\overline{L}_{\text{pruned}})}{\sqrt{np}(abp)}=f(np)$$
(23)
for $`n1`$, where $`a=0.872\pm 0.002`$, $`b=0.105\pm 0.005`$ and $`f(np)1`$ for large $`np`$. The shift factors indicate that the PTSP strategy can be no more than $`\frac{0.872}{0.767}1=14\left(\pm 1\right)`$ worse than the re-optimization strategy.

The master curve for the a priori tour length is shown in fig. 4. Our scaling arguments predict that the shift factors $`\beta _{\text{a priori}}(p)`$ should tend towards a constant for $`p0`$. However, data are fit very well by the relation

$$\beta _{\text{a priori}}(p)=\frac{1}{1.250.82ln(p)}$$
(24)
which would tend to zero as $`p0`$ in conflict with our scaling arguments. To resolve this dilemma we need to probe very small $`p`$.

## V The Limiting Case $`p0`$

We are interested in finding whether $`\beta _{\text{a priori}}(p)`$ tends towards a constant as $`p0`$. To do this using the above approach is difficult, since we need a large number of cities to produce reliable data for this regime. Extraction of this behaviour can however be achieved by comparing simulations for different values of $`n`$, but fixed $`np`$. We accomplish this by insisting that each instance has 4 cities on the pruned tour. 4 city tours are chosen since they are the smallest for which it matters in which order the cities are visited. This can be viewed as an efficient way to simulate (approximately) the PTSP strategy with $`p=\frac{4}{n}`$.

Since we are considering the PTSP at fixed $`np`$, if $`\beta _{\text{a priori}}(p)`$ tends towards a (non-zero) constant as $`p0`$ then we expect $`E\left(L_{\text{a priori}}^{\text{4 city}}\right)/n`$ to tend towards a constant as $`p0`$. Simulations in this regime were performed for $`N=12210`$, with $`100`$ different random city configurations used for $`N<30`$, $`20`$ configurations for $`N90`$ and $`10`$ configurations for $`N120`$. Figure 5 shows a linear-log plot of $`\frac{n}{2E\left(L_{\text{a priori}}^{\text{4 city}}\right)}`$ against $`\mathrm{ln}(n/4)\mathrm{ln}(1/p)`$. For small $`n`$ these results reasonably match the direct measurements of $`\beta _{\text{a priori}}(p)`$, shown for comparison. However, for surprisingly large $`n100`$ which is beyond the range of our $`\beta _{\text{a priori}}(p)`$ data, our earlier proposal of scaling behaviour is vindicated by $`E\left(L_{\text{a priori}}^{\text{4 city}}\right)/n`$ approaching a constant value. In summary we have

$$E(L_{\text{a priori}})=\sqrt{\frac{n}{p}}\beta _{\text{a priori}}(p)$$
(25)
where

$$\beta _{\text{a priori}}(p)\{\begin{array}{cc}=\frac{1}{1.250.82ln(p)}\hfill & p>0.03\hfill \\ =\beta _0\hfill & p<0.03.\hfill \end{array}$$
(26)
## VI Notes on Algorithm Implementation

We applied stochastic annealing to the PTSP using a combination of the 2-opt and 1-shift move-setslin1  established for the TSP. Both move-sets work similarly to that which would be expected for the deterministic case. The expected pruned tour length change for the move was estimated by averaging the change in the tour length for a number of instances. For a given instance it is not necessary to decide whether every city is present, but only the set of cities closest to the move which determine the change in the pruned tour length (see figure 6). For the PTSP, the location of the nearest cities on the pruned tour to the move is determined from a simple Poisson distribution.

When using stochastic optimization, the only variable over which we have control is the sample size (the number of instances) $`r`$, whereas the effective temperature $`\frac{\sigma }{\sqrt{r}}`$ also entails the standard deviation $`\sigma `$ of the pruned length change over instances. As shown in figure 7, annealing by controlling $`r`$ alone exhibits a relatively sharp transition in the expected pruned tour length. The rapid transition appears to ‘freeze in’ limitations in the tours found (analogous to defects in a physical low temperature phase). By comparison we obtain a much smoother change when $`\frac{\sigma }{\sqrt{r}}`$ is controlled.

The sharpness of the transition under control by $`r`$ is caused by the fact that $`\sigma `$ may vary from move to move, and is on average lower when the expected pruned tour length is less. The jump in the pruned tour length is accompanied by a jump in $`\sigma `$ and hence the temperature. We suggest that quite generally controlling $`\frac{\sigma }{\sqrt{r}}`$ gives a better cooling schedule than focussing on $`r`$ alone.

## VII Conclusion

We have shown that earlier incompatible ideas about the form of PTSP tours especially at small $`p`$bertsimas1 ; bertsimas4 ; bartholdi  are resolved by a new crossover scaling interpretation. The crossover scale corresponds to a group of cities such that of order one will typically have to be visited; below this scale the (optimsed) a priori PTSP tours resemble a local sort whereas they are TSP-like on scales larger than the crossover. Our computational results for the pruned tour length are summarised by eq. 23 and clearly support the crossover scaling.

Computationally the a priori tour length is more subtle than the pruned tour length, although it does ultimately conform to expectations from crossover scaling. We introduced 4-city tours to probe the behaviour of a priori tour length down to very small $`p`$. As summarised by eq. 25, we find a wide pre-asymptotic regime until recovering the expected crossover scaling only for $`p<0.03`$. Understanding these anomalies in the a priori tour length, and confirming them analytically, is left as a future challenge.

We have shown stochastic annealing to be a robust and effective stochastic optimization technique, taking the PTSP as a representative difficult stochastic optimization problem. In this case it enabled us to obtain representative results out to unprecedented problem sizes, which in turn supported a whole new view of how the tours behave. Of relevance to wider applications of stochastic optimmisation, we have seen that smoother annealing can be obtained by directly controlling the effective temperature$`\frac{\sigma }{\sqrt{r}}`$fink1  rather than simply the bare depth of sampling $`r`$ alone.

NEB would like to thank BP Amoco & EPSRC for the support of a CASE award during this research.

## The Length of a TSP Tour for small numbers of cities

Numerical estimates of the length of a TSP tour for $`n10`$ are given below


text:
# An Extinction Map and Color Magnitude Diagram for the Globular Cluster NGC 3201

## 1. Introduction

We are currently undertaking a survey of approximately 15 Galactic globular clusters (GCs) with the aim of identifying eclipsing binary stars around or below the turn-off-point by means of detecting brightness variations. These systems can provide masses for Population II main sequence stars. High-quality photometry is essential for obtaining reliable values for surface temperatures, luminosities, ages, etc. Precise extinction determinations are a critical part of obtaining these results from our photometry.

NGC 3201, by chance one of the first GCs we analyzed, is a low-latitude cluster located at $`\alpha _{2000}=10^h36^m36.8^s`$ and $`\delta _{2000}=46^{}24^{^{}}40^{^{\prime \prime }}`$ or $`l=277.2^{}`$ and $`b=8.6^{}`$ (Harris 1996) with a retrograde orbit around the center of the Milky Way. Its nearby location ($`d`$ = 5.2 kpc) and low concentration ($`c`$ = 1.31) (Harris 1996) make it an attractive target for photometric studies. Recent such studies include Cacciari (1984), Brewer et al. (1993), Côté et al. (1994), Covino & Ortolani (1997) (hereafter C&O), and Gonzalez & Wallerstein (1998) (hereafter GW).

Due to its low-latitude position, the effects of differential reddening across the face of NGC 3201 are quite substantial. The existence of variable extinction was noted in practically all earlier studies of this cluster. Cacciari (1984), for instance, finds an irregular reddening distribution with a variation of $`\mathrm{\Delta }E_{BV}`$ = 0.03 mag in addition to a mean value of $`E_{BV}`$ = 0.21. GW report a range in $`E_{BV}`$ of as much as 0.1 mag across NGC 3201. Using spectroscopic data for 18 giants in NGC 3201, GW correct for the color spread by modelling the reddening in $`E_{BV}`$ as a plane in $`\alpha `$ and $`\delta `$.

Recently published dust infrared emission maps by Schlegel, Finkbeiner, & Davis (1998) (SFD hereafter) seem to indicate, however, that the dust distribution in the region of NGC 3201 is too clumpy to be fitted by a linear function in position, even on the scale of arcminutes (see Fig. 6).

In this work, we attempt to remove the differential reddening across the cluster by using our high-quality $`VI`$ photometry. Our procedure is, in principle, similar to the one used by Piotto et al. (1999). It aims to find the average $`\mathrm{\Delta }E_{VI}`$ for subregions of the cluster field of approximately 1 arcmin2 in size with respect to a fiducial region in the cluster where little or no differential reddening is apparent and where the overall $`E_{VI}`$ is small compared to the rest of the cluster. The pixel size (resolution) of this extinction map is approximately one fourth of the size (in area) of the SFD maps. Since the SFD maps are based on infrared emission from dust, our map will provide an independent check of whether, if at all, they tend to overestimate reddening, as was suggested by, e.g., Arce & Goodman (1999) for regions with $`A_V>0.5`$ mag.

Our results indicate variations in $`E_{VI}`$ of up to $``$ 0.2 mag. Differential reddening this strong can wreak havoc with photometric and spectroscopic studies of cluster stars. The inherent uncertainties in the various parameters of the CMD (e.g., magnitude and color of the main-sequence turnoff (MSTO), horizontal branch features, etc.) are greatly amplified. Moreover, one may “detect” age or surface-temperature gradients where there are none in reality. To give an example of the magnitude of the effect, using the recently published color-temperature relations by Houdashelt, Bell, & Sweigart (2000), the effective temperature of a solar-metallicity, main sequence star with $`VI0.8`$ would vary by approximately 600$`K`$ for a differential reddening effect of $`E_{VI}0.2`$ mag, which, in turn, could lead to errors in the metallicity determination.

Using our internal dereddening technique, we obtain a high-quality, deep-photometry CMD of NGC 3201 comprised of approximately ninety 600s exposures as well as some shorter ones, all using the $`V`$ and $`I`$ bands.

Our observations and data reductions, as well as the details about our internal dereddening method, are documented in Section 2. Section 3 contains our results concerning the reddening map and the CMD of NGC 3201. Finally, we discuss the determination of the reddening zeropoint of our extinction map in Section 4 and give a brief summary of our work in Section 5.

## 2. Observations, Data Reduction, and Internal Dereddening

### 2.1. Observations and Basic Data Reductions

The NGC 3201 observations were obtained during the nights of April 26 through May 6, 1998, at the Las-Campanas-Observatory (LCO) 1m Swope Telescope, using Johnson-Cousins $`VI`$ filters and a SITe 1 $`2048^2`$ CCD with a field-of-view of 23.5 arcmin on a side. Table 1 gives the number of epochs we observed for different exposure times which were chosen in order to cover a larger magnitude range in the CMD.

The initial processing of the raw CCD images was done with the routines in the IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF. CCDPROC package. For each night, 10 bias frames were combined for the bias subtraction. The $`V`$ band flats were produced by combining between 4 and 6 twilight flat images per night. All $`I`$ band data were first flattened using an image comprised of 15 individual $`I`$ band domeflats, and then divided by a normalized, dome-flattened dark-sky-flat which itself was created by median-averaging approximately 40 individual different blank sky fields. Finally, the $`I`$ band images were corrected for fringing by subtracting a fringe image which was created by subtracting the mean pixel value of the dark-sky-flat from the dark-sky-flat image itself.

The processed data were reduced using DoPHOT (Schechter, Mateo, and Saha 1993) in “fixed-position mode” where the positions of stars are fixed (after correcting for global frame-to-frame shifts as well as small distortions) to the positions measured on a deep-photometry template image obtained by co-adding the $``$ 15 best-seeing frames.

The photometric results for every star were averaged over all frames with the same exposure time to obtain the final magnitude for the star under investigation. For the 600s exposures, our only requirement was that a star appear in more than 75% of the epochs. The 60s images were taken to complete the CMD in the brighter regions; we only took magnitudes from these exposures of stars which were saturated in the 600s exposures. The same procedure was followed for the 10s frames. Aperture corrections were applied by calculating one constant value for the whole chip after we found little dependence upon position on the CCD. The data from the shorter exposures were shifted to the photometric system of the 600s exposures using non-saturated stars in common.

During our consecutive two photometric nights (May 2 & 3, 1998), we obtained a total of more than 140 observations of various Landolt (1992) standard stars which covered a range of $`0.2<VI<2.2`$ in color as well as $`1.15<X<1.55`$ in airmass. Using the IRAF PHOTCAL package, we applied one single standard star solution for data from both nights of the form

$$V=v+a_0+b_0X_v+c_0(vi)$$
(1)
$$VI=a_1+b_1X_i+c_1(vi),$$
(2)
where the $`a_j`$, $`b_j`$, and $`c_j`$ are the fitted constants, $`X_{filter}`$ is the airmass of the exposure taken with the respective filter, the lowercase magnitudes are instrumental, and the uppercase ones are the known magnitudes of the standard stars. The root mean square error was 0.007 mag and 0.016 mag for $`V`$ and $`VI`$, respectively

Astrometry was performed by identifying 59 USNO reference stars (Monet et al. 1996) in NGC 3201 and using the IRAF IMAGES package for the coordinate transformation. A linear fit in $`x`$ and $`y`$ produced errors around 0.2 arcsec, consistent with the USNO precision. The tangent point of the transformation was $`(x,y)=(962.7,934.2)`$ and corresponds to $`\alpha _{2000}=10^h17^m31.2^s`$ and $`\delta _{2000}=46^{}24^{^{}}11^{^{\prime \prime }}`$. The rotation of $`x`$ and $`y`$ with respect to $`\alpha `$ and $`\delta `$ was $`91.23^{}`$, and the pixel scale came out to be 0.695 pix/arcsec. These results enabled us to create the reference grid in Fig. 3.

### 2.2. Mapping the Differential Reddening

At this point in the data reduction we were in the possession of a fully calibrated CMD whose main sequence appeared exceptionally broad (see Fig. 8). After ruling out a number of possible instrumental causes for this effect, and after we noted a clear dependence of the appearance of the main sequence upon the positions of stars on the CCD (illustrated in Fig. 1), we concluded that the underlying cause for the broadness of our main sequence was differential extinction across the field of view. This seemed especially sensible given the low-latitude location of NGC 3201 and motivated us to create a reddening map since our binary star data were not useful without this correction. In order to create such a map for NGC 3201, we developed the steps described below (see also Fig. 2):

1. A fiducial region was chosen in which

* very little or no differential reddening occurred, i.e., where the main sequence appeared narrow;
    * the overall $`E_{VI}`$ was very small with respect to the rest of the field of view;
    * there were enough stars in the field to obtain a statistically reliable fit to their positions on the CMD.

The coordinates of this fiducial region are approximately $`10^h16^m26^s<\alpha _{2000}<10^h16^m53^s`$ and $`46^{}31^{^{}}21^{^{\prime \prime }}<\delta _{2000}<46^{}26^{^{}}44^{^{\prime \prime }}`$ (see Figures 3 through 7).
2. A high-order polynomial was fit to the main sequence stars with $`18<V<21`$ and $`0.7<VI<1.4`$ in the fiducial region, using algorithms supplied by Robbie Dohm-Palmer (2000, private communication) which enabled us to separate the main sequence stars from background stars in the CMD, and by Press et al. (1992), which provided the least-squares fit to the data (see Fig. 2).
3. We divided the field of view into subregions of different sizes in order to have a sufficient number of stars in each region. Towards the outer parts of the field, the stellar density decreases. As a result, the sizes of the subregions increase. The inner subregions are all of the same size (100$`\times `$100 pixel). See Fig. 3.
4. For each of these regions, every star falling between $`17.9<V<21.1`$ and $`0.65<VI<1.45`$ was incrementally moved along the reddening vector defined by the relations in Cardelli, Clayton, & Mathis (1989) until it intersected the fit (see Fig. 2).
5. The statistical biweight (see Beers, Flynn, & Gebhardt 1990) of all these incremental shifts of the stars in a given subregion was then calculated, outliers (0.9 $`\sigma `$) were removed, and it was recalculated. This provided the value of the total shift along the reddening vector for each star. Since the slope of the reddening vector for a standard extinction law is known (1.919; see Cardelli, Clayton, & Mathis (1989)), the differential $`E_{VI}`$ for each star corresponds to the $`(VI)`$-component of the vector described above (Fig. 2).
6. For six of the 100$`\times `$100 pixel subregions as well as one 200$`\times `$300 pixel subregion on the extinction map shown in Fig. 3, too few stars were present to calculate the average reddening. For these cases, we calculated $`E_{VI}`$ by simply averaging the reddening values of all the neighboring subregions222In this process, every 100$`\times `$100 subregion was treated individually, even if it was a part of a larger one. That is, for the 200$`\times `$300 region just south of the north-east corner of the field of view in which the number of stars was too low to obtain a reddening value, the above procedure was applied six times..
7. Our error analysis for the shifts described above involved the following:

* We created an error ellipse for every star defined by the values of the associated errors (as returned by DoPHOT) in color and magnitude.
    * Since the color and magnitude errors are correlated, the error ellipse is tilted, with the tilt angle and the lengths of the semimajor and semiminor axes functions of $`\sigma _V`$ and $`\sigma _{VI}`$.
    * One average tilt angle was obtained for every subregion under investigation by averaging the results for the individual stars. For every subregion, this angle was approximately $`35^{}\pm 4^{}`$. The semimajor and semiminor axes for every subregion were obtained in the same way, so that every star in a given subregion has the same error ellipse associated with it in the center.
    * This error ellipse was shifted along with its respective star during the dereddening process described above in items 4 and 5. The “point of first contact” (pfc), i.e, the point at which the error ellipse first touches the fit through the data in the fiducial region, as well as the “point of last contact” (plc), represent the respective one-sigma-deviation points (see Fig. 2).
    * These two contact points are not necessarily symmetric about the reddening value of the star (center of error ellipse), but in order to increase the readability of this publication, we calculated the mean pfc and plc for a given subregion, and then their average distance from the center of the ellipse, so that $`\sigma _{ellipse}=0.5(pfcplc)`$.
    * The final error estimate for each subregion ($`\sigma _{E_{VI}}`$) given as the lower number in each of the pixels in Fig. 4 is then obtained by adding in quadrature $`\sigma _{ellipse}`$ and the error in the mean of the shifts of all the stars in the subregion under investigation.

## 3. Results

### 3.1. Extinction Map and Comparison with SFD Map

Fig. 3 summarizes the results of our internal dereddening procedure. Darker regions correspond to higher reddening with respect to the average reddening taking place in the fiducial region described in Section 2.2, item 1. For reference purposes, we included the approximate locations of this fiducial region and the core radius of NGC 3201, as well as a coordinate grid. The average $`E_{VI}`$ relative to the fiducial region is 86 mmag with a standard deviation of 61 mmag. The individual subregions’ $`E_{VI}`$ values in millimagnitudes (relative to the average reddening in the fiducial region) are shown as the top number in each of the pixels in the grid of Fig. 4. The bottom number in each pixel corresponds to the error in the mean for the corresponding $`E_{VI}`$ value. In addition, we included an image of NGC 3201 (Fig. 5) in order to indicate the locations of the fiducial region and the core radius of the globular cluster.

At this point, these reddening values are all relative to our fiducial region which itself is suffering some mean interstellar extinction. The procedure described cannot determine the absolute reddening zeropoint. Thus, we had to introduce additional information, derived from previous studies, in order to calculate this zeropoint. The discussion of this process is given in Section 4.1.

As a first step, we subtracted our extinction map from the SFD map of the region of NGC 3201 which we observed. The SFD map of the cluster region is shown in Fig. 6 with the same orientation and size as Fig. 3. Once again, dark corresponds to higher extinction. We added the location of the fiducial region for reference. The average $`E_{VI}`$ of the SFD map is 487 mmag with a standard deviation of 59 mmag.

A first comparison indicates that, while the larger-scale structure is the same for both extinction maps (in particular the obscuration ridge extending from the top left part of the image to the right center part of it), our map shows smaller-scale features which are not present in the SFD maps. These features are visible in the difference image in Fig. 7 whose orientation and size is once again identical to the other Figures. Darker regions correspond to places where our data indicate the presence of dust which does not show up in the SFD map; lighter regions show good agreement between the two maps. The average $`E_{VI}`$ of this difference map is 401 mmag with a standard deviation of 46 mmag (significantly smaller than the standard deviation of both our extinction map and the SFD map). This lower standard deviation is an indication that both our map and the SFD map of the region of NGC 3201 are tracing many of the same features in the foreground extinction along the line of sight to the globular cluster.

The smaller-scale features present in only our map but absent in the SFD map are certainly worth further attention. At this point, we are not entirely confident that we can give a definite reason for the discrepancies, but some potential reasons could be the following:

* The smaller-scale features which are not in the SFD maps might have been smoothed out by the lower resolution of the IR observations used to create the maps.
* The temperature of the dust producing the smaller-scale features was outside the detection range of the instrument used to collect the data for the SFD maps. This way, they would still show up in our analysis.
* The dust we are seeing is actually part of the globular cluster itself333The phenomenon of dust inside globular clusters was discussed by, e.g., Forte et al. (1992). and was therefore not picked up in the SFD maps.

### 3.2. Photometry and CMD

Fig. 8 contains the raw CMD (no dereddening applied) for NGC 3201. To illustrate the effects of our dereddening method, we show, in Fig. 9, the CMD after applying the extinction map (differential reddening correction) to the data. The datapoints are now all shifted to the CMD-location of the stars in the fiducial region, that is, no absolute reddening zeropoint has been applied. The improvement is quite stunning. The width of the main sequence of Fig. 9 has, by applying the differential reddening map, decreased to a fraction of its former value in Fig. 8. The subgiant branch which was hardly visible at all in the raw CMD is now well traced out by the datapoints. Even the horizontal branch, though sparsely populated in both CMDs due to saturation effects of the bright stars in the field, is narrower in the differentially dereddened CMD.

## 4. The Reddening Zeropoint of our Extinction Map for NGC 3201

### 4.1. Comparing the MSTO Color with Previous Studies

In order to make our reddening map (Fig. 4) a useful tool, we need to determine the reddening zeropoint to add to the $`E_{VI}`$ values in the grid. As mentioned above, this could not be done using our data alone. In order to be able to use results from previous studies for the calculation of this zeropoint, an agreement between our photometry results and the corresponding literature values is necessary. When comparing MSTO color, however, one needs to consider that different parts of the field of view are reddened by different amounts, as we showed in our extinction map (Fig. 3).

We estimated our main-sequence turnoff $`V_{TO}`$ magnitude to be 18.2 and the MSTO color $`(VI)_{TO}`$ to be around 0.88 mag; both were determined by eye from the raw photometry CMD (Fig. 8).

These values were compared to the results of Rosenberg et al. (2000) (hereafter RB) and C&O. RB find $`V_{TO}18.2`$ and $`(VI)_{TO}0.92`$; both values are again determined by eye. C&O report $`V_{TO}18.15`$ and $`(BV)_{TO}0.65`$ (tabulated). Using their reddening estimate of $`E_{BV}0.22`$, the relations in Cardelli, Clayton, & Mathis (1989), and a set of isochrones provided by Don VandenBerg for $`[Fe/H]=1.41`$ (VandenBerg 2000, private communication, based on evolutionary models by VandenBerg et al. 2000; hereafter VDB), the C&O value for $`(VI)_{TO}`$ can be calculated to be approximately 0.93.

One discrepancy between our results and these literature values is therefore that our value for $`(VI)_{TO}`$ is approximately 0.05 mag bluer than what is quoted by the authors mentioned above. The existence of variable reddening may be a likely explanation for this discrepancy, especially when one considers which part of the cluster was observed.

RB concentrate on the southern and eastern parts of the cluster for their longer exposures and the center of NGC 3201 for their shorter exposures. All of these three regions are, according to our reddening map (Fig. 3), regions of higher extinction with respect to our fiducial region.

C&O observe regions just to the north and west of the cluster center and the cluster center itself. While the north and west regions suffer only low reddening (Fig. 3), the area within $``$ 7 arcmin to the north and west of the cluster center is obscured by a smaller feature not visible in the SFD map (see Fig. 6). At least 5 of 7 C&O fields are very likely affected by this dense obscuration feature.

While a quantitative (star-by-star) magnitude comparison between our data and the data obtained by RB and C&O is not possible, we conclude that their estimates for the MSTO color could very well be influenced by the fact that the regions they observed were subject to extinction somewhat higher than the average value over our field of view.

We note that these comparisons focus on the MSTO colors; the vertical nature of the MSTO makes it difficult to compare the $`V`$ mag of the MSTO at sufficiently high precision.

### 4.2. Basic Isochrone Fitting and Determination of Reddening Zeropoint

We examined three different possibilities of finding the $`E_{VI}`$ zeropoint to add to our extinction map. All three methods are outlined below. For each case, it is worth pointing out that a direct comparison with literature values should be taken with caution since usually, only a single numerical value for the reddening is given, whereas our result is a reddening matrix of which every single element needs to be added to a zeropoint offset.

#### 4.2.1 Using Isochrone Fits

Simultaneously fitting age, distance, and reddening to the aforementioned set of VDB isochrones produced the fit shown in Fig. 10. It should be noted that since the point of this publication is outlining our dereddening method as well as creating an extinction map for NGC 3201, we only performed some basic isochrone fitting here. For an estimate of $`[Fe/H]`$, we adopted the isochrone with metallicity closest to the GW value of -1.42.

For the fit of Fig. 10, $`[Fe/H]=1.41`$, $`d4.5`$ kpc, age = 18 Gyrs, and the $`E_{VI}`$ zeropoint to be added to the values in Fig. 4 is 0.15 mag. The age is clearly on the high end of things, but the 18 Gyr isochrones produced by far the best fits, even when $`[Fe/H]`$ and the distance were slightly modified. The value for the distance is slightly below the one quoted by Harris (1996), but fairly close to the two values derived using slightly different methods by C&O.

When adding the mean $`E_{VI}`$ of our reddening map (86 mmag) to the zeropoint, the average $`E_{VI}`$ for the cluster is approximately 0.24, which is slightly below the estimates of, e.g., Cacciari (1984) whose $`E_{VI}0.32`$, and of Harris (1996) and RB who quote the same value.

In order to justify our choice of isochrone parameters and to eliminate the possibility of multiple factors conspiring against us in the attempt to find a good isochrone fit to the data, we illustrate in Fig. 11 what the CMD would look like if we were to strictly adopt literature values from Harris (1996) ($`[Fe/H]=1.58`$, $`E_{VI}=0.325`$, d = 5.2 kpc). Plotted are a pair of isochrones with respective ages 12 and 14 Gyrs on top of the data444In order to have a valid basis for comparison, we needed to apply an $`E_{VI}`$ zeropoint of $`0.3250.086=0.239`$ to the differentially dereddened data, where 0.325 is the Harris (1996) $`E_{VI}`$ value and 0.086 is the mean $`E_{VI}`$ value of our extinction map.. The agreement between data and fit is not as good as in Fig. 10.

#### 4.2.2 Using SFD Maps

As stated in Section 3.1, the mean $`E_{VI}`$ of the difference map between the SFD and our data is $`401\pm 46`$ mmag. The average $`E_{VI}`$ value for this difference map corresponds to the $`E_{VI}`$ zeropoint to add to our data to obtain absolute reddening values.

Calculating the total reddening in this fashion gives $`E_{VI}0.49`$ which is significantly higher than the literature values. Furthermore, using the isochrones mentioned in the previous subsection, we were not able to produce a fit to the data due to the fact that $`E_{VI}`$ was too high.

#### 4.2.3 Using RRLyrae Stars

Another method of attempting to find a reddening zeropoint is to simply adopt one from the literature and model the thus dereddened data with a set of isochrones. Here, we took Cacciari’s (1984) value of $`E_{VI}0.32`$ which she calculated using RRLyrae stars. Since this corresponds to the average reddening for the whole cluster, we had to subtract the mean $`E_{VI}`$ of her RRLyrae stars (in our field of view) with respect to our fiducial region. This value turned out to be 79 mmag.

We failed, however, to find a good isochrone fit to the thus dereddened data, as is shown in Fig. 12. In order to get the MSTO points of the data and the isochrone to agree, we had to lower $`[Fe/H]`$ to -2.01, a value quite a bit below literature values such as -1.41 (GW) or -1.58 (Harris 1996). The distance modulus here is 13.4, and the age 14 Gyrs. As one can easily see, the fit fails to hit the subgiant branch due to the “young” isochrone age of 14 Gyrs. Older age isochrones could not be shifted far enough toward the dereddened MSTO location to produce a fit, even when the metallicity was lowered to the lowest available value of -2.31.

#### 4.2.4 Comments on the Adopted Reddening Zeropoint

Given that we were not able to produce a fit to the data using the VDB isochrones when we calculated the reddening zeropoint using either the SFD maps or the Cacciari (1984) RRLyrae stars, we conclude that our best estimate for the $`E_{VI}`$ zeropoint is 0.15 mag, giving an average $`E_{VI}0.24`$ for the cluster as a whole. This value falls approximately $`1.5\sigma `$ below the estimate by Cacciari (1984) and the value tabulated by Harris (1996).

We could force our reddening zeropoint to agree with these past results by applying a systematic shift to our $`I`$ magnitudes of 0.08 to brighter values which would give a distance to NGC 3201 of 4.8 kpc. We have, however, no reason to adopt such a shift. The transformations to standard magnitudes as well as the application of the aperture corrections (see Section 2.1) underwent repeated thorough examinations and appear to be correct. Since no other indication throughout our data reduction and analysis suggests a systematic error in the $`I`$ band magnitudes, either, we simply cannot justify adjusting our $`I`$ magnitudes in order to compensate for the $`1.5\sigma `$ effect described above.

## 5. Summary and Concluding Remarks

In the process of finding eclipsing binary stars in NGC 3201, we noticed the existence of variable reddening of up to 0.2 mag in $`E_{VI}`$ on a scale of arcminutes. Using our internal dereddening method outlined in Section 2.2, we created an extinction map which is shown in Figures 3 and 4. Applying the map to our raw data (Fig. 8) significantly improved the appearance of the CMD (Fig. 9).

Comparison between our extinction map and the SFD map of the same region (Fig. 6) showed that the same larger-scale features exist in both maps. Our map displays some additional, smaller-scale features which are absent in the SFD maps (see Fig. 7).

The $`E_{VI}`$ zeropoint which needs to be added to the numbers in Fig. 4 in order to get absolute $`E_{VI}`$ values is 0.15. This value is below literature results (Cacciari 1984, Harris 1996) by approximately $`1.5\sigma `$, but produced by far the best VDB isochrone fit to the data. The zeropoint determined with the help of the SFD maps gives $`E_{VI}0.49`$ as the average value across NGC 3201 which is higher than literature values which supports the statement by Arce & Goodman (1999) that the SFD maps overestimate the reddening in regions of high extinction.

The results from this work will be essential in our binary star research where high-quality photometry of every binary system is necessary for distance determinations. A vital condition to obtaining these measurements is of course the knowledge of the exact extinction the star under investigation is suffering.

Furthermore, studies like this may be useful in determining properties of the ISM itself such as examining a possible dependence of $`R_V`$ upon position in the field of view or giving insight into the distribution and properties of the dust along the line of sight.

This research was funded in part by NSF grants AST96-19632 and AST98-20608. We would like to sincerely thank Alex Athey, Kristin Chiboucas, Robbie Dohm-Palmer, Ivan King, and Don VandenBerg, as well as the anonymous referee, for their help, advice, and pointing out errors during the various stages of the creation of this publication.


text:
# 1 Linear Algebra

## 1 Linear Algebra

In this paper, the bold-faced characters are reserved for denoting vectors; the dimension of the vectors is $`n3`$, except when the context suggests otherwise. For example, $`𝐟`$ is a vector with coördinates $`f_1,\mathrm{},f_n`$, while $`\mathrm{𝟎}`$ and $`\mathrm{𝟏}`$ are the all-zero and all-one vectors, respectively.

Consider the system of linear equations

$$(1)A𝐱=𝐟$$

where $`𝐱`$ is the vector of the unknowns, $`𝐟`$ is a fixed vector, and $`A`$ is an $`n\times n`$ matrix such that:

1. $`A`$ is symmetric;
2. the off-diagonal entries of $`A`$ are non-positive;
3. $`A`$ is irreducible, which is to say that the simple graph whose adjacency matrix has the same off-diagonal zero pattern as $`A`$ is connected (throughout the paper, we refer to this graph as “the underlying graph $`G`$”);
4. $`A\mathrm{𝟏}=\mathrm{𝟎}`$.

Applied to a matrix with properties (i)–(iii), the Perron-Frobenius Theorem states that its least eigenvalue has multiplicity one, and the corresponding eigenvector can be chosen to have all its coördinates positive. Hence, by (iv), $`A`$ is positive semidefinite and of rank $`n1`$. Consequently, (1) has no solution if $`𝐟\mathrm{𝟏}0`$ and a solution unique up to a shift by a multiple of $`\mathrm{𝟏}`$ if $`𝐟\mathrm{𝟏}=0`$. For $`i\{1,\mathrm{},n\}`$, set $`𝐟_i:=(f_1,\mathrm{},f_i𝐟\mathrm{𝟏},\mathrm{},f_n)`$ (so that $`𝐟_i\mathrm{𝟏}=0`$) and let $`𝐱_i=(x_{i1},\mathrm{},x_{in})`$ denote the \[unique\] solution to the system of linear equations

$$(2)A𝐱=𝐟_i$$

such that $`x_{ii}=0`$. Our goal is to demonstrate the identity

$$(3)x_{12}+x_{23}+x_{31}=x_{21}+x_{32}+x_{13}.$$

Suppose first that $`n>3`$ and let us apply \[symmetric\] Gaussian elimination to (2), to eliminate the non-zero off-diagonal entries of $`A`$ in the last column and row. The reader can easily see that the first $`n1`$ unknowns are not affected by such an elimination, while the upper left $`(n1)\times (n1)`$ corner, $`A^{}`$, of the matrix obtained from $`A`$, and the $`(n1)`$-dimensional vector, $`𝐟_i^{}`$, whose entries are the same as the first $`n1`$ entries of the vector obtained from $`𝐟_i`$, have the same properties as $`A`$ and $`𝐟_i`$, respectively. In other words, the system of $`n1`$ linear equations

$$(4)A^{}𝐱^{}=𝐟_i^{}$$

with $`n1`$ unknowns $`x_1,\mathrm{},x_{n1}`$ can be considered in lieu of (2), and if $`𝐱_i=(x_{i1},\mathrm{},x_{in})`$ is a solution to (2) then $`𝐱_i^{}:=(x_{i1},\mathrm{},x_{i,n1})`$ is a solution to (4). Hence, in order to establish (3), we may assume that $`n=3`$.

Still, to verify (3) for $`n=3`$ directly is a rather tedious task. After all, one would have to solve three systems of three equations each in general. Instead, we prefer the following detour. In the next section, we will see that (3) is trivially true if the underlying graph $`G`$ is the star $`K_{1,n1}`$. If $`n=3`$ then $`G`$ is either the star $`K_{1,2}`$ or the triangle $`K_3`$. In the latter case, the matrix $`A`$ has the form

$$\left[\begin{array}{ccc}\alpha _2+\alpha _3& \alpha _3& \alpha _2\\ \alpha _3& \alpha _1+\alpha _3& \alpha _1\\ \alpha _2& \alpha _1& \alpha _1+\alpha _2\end{array}\right]$$

where $`\alpha _1,\alpha _2,`$ and $`\alpha _3`$ are some positive numbers. Set $`c:=\alpha _1\alpha _2+\alpha _2\alpha _3+\alpha _3\alpha _1`$ and $`\beta _i:=c/\alpha _i`$, $`i\{1,2,3\}`$. Then, the matrix

$$B:=\left[\begin{array}{cccc}\beta _1& 0& 0& \beta _1\\ 0& \beta _2& 0& \beta _2\\ 0& 0& \beta _3& \beta _3\\ \beta _1& \beta _2& \beta _3& \beta _1+\beta _2+\beta _3\end{array}\right]$$

has properties (i)–(iv) and it is not difficult to verify that $`A`$ is the upper left $`3\times 3`$ corner of the matrix obtained from $`B`$ by symmetric Gaussian elimination. Hence, (3) holds for $`A`$ if it holds for $`B`$. Which it does, since the underlying graph for $`B`$ is the star $`K_{1,3}`$.

## 2 Mechanical Model

The systems (1) and (2) have the following mechanical interpretation: $`n`$ masses labelled $`1`$ through $`n`$ are submerged in water, so that each mass $`j`$ is subject to a (positive or negative or zero) weight, $`f_j`$. In addition, every two masses $`j`$ and $`k`$ are connected by a rubber band of resistance $`a_{jk}=a_{kj}0`$ where $`a_{jk}`$ is the $`jk`$ entry of $`A`$. (A rubber band is of resistance at most $`a0`$ if and only if a force of magnitude $`a\mathrm{}`$ is sufficient to stretch it to length $`\mathrm{}`$.) If mass $`i`$ is nailed to the origin, such a system of masses has a unique equilibrium. Let $`x_{ij}`$ denote the altitude of mass $`j`$ in this equilibrium. (So that $`x_{ii}=0`$.) The fact that the system is in an equilibrium is expressed by equating the sum of the forces acting on every mass $`j`$ to zero. In an equilibrium, these forces are the weight $`f_j`$, the reaction $`𝐟\mathrm{𝟏}=(f_1+\mathrm{}+f_n)`$ of the nail if $`j=i`$, and the resistance forces $`a_{1j}(x_{i1}x_{ij}),\mathrm{},a_{nj}(x_{in}x_{ij})`$. In other words, the vector $`(x_{i1},\mathrm{},x_{in})`$ is the unique solution $`𝐱_i`$ to the system (2) such that $`x_{ii}=0`$.

Using this physical interpretation of (2), we see that (3) is trivially true if the underlying graph $`G`$ is the star $`K_{1,n1}`$ with say center $`k`$. This is because $`x_{ij}=x_{ik}+x_{kj}`$ in this case. This completes the proof of (3).

Remark 1:  In the mechanical model, we could replace the weights $`f_i`$ by arbitrary $`d`$-dimensional forces and, correspondingly, the altitudes $`x_{ij}`$ by the $`d`$-dimensional position vectors. Then, the “$`d`$-dimensional generalization” of (3) would be proved coördinate-wise.

Remark 2:  The identity $`x_{i_1i_2}+x_{i_2i_3}+\mathrm{}+x_{i_ki_1}=x_{i_2i_1}+x_{i_3i_2}+\mathrm{}+x_{i_1i_k}`$ can be shown for any sequence $`i_1,\mathrm{},i_k`$ of indices using (3); we leave this to the reader.

## 3 Random Walks

Let $`G`$ be a connected graph on $`n3`$ vertices. To obtain a random walk out of vertex $`j`$ in $`G`$, one recursively builds a sequence $`i_0i_1i_2\mathrm{}`$ of vertices, in which $`i_0=j`$ and, for $`t0`$, the vertex $`i_{t+1}`$ is chosen among the neighbors of $`i_t`$ with the uniform probability. The hitting time $`H(j,i)`$ is the expectation of the least number $`t`$ such that $`i_t=i`$. (So that $`H(i,i)=0`$.) If $`ij`$ then, clearly,

$$(5)H(j,i)=1+\frac{1}{𝚍eg(j)}\underset{\{k:jkG\}}{}H(k,i).$$

Multiplying both sides of (5) by $`𝚍eg(j)`$, we see that the vector $`𝐱_i:=(H(1,i),\mathrm{},H(n,i))`$ is a solution to the system of linear equations

$$(6)A𝐱=𝐠_i$$

where $`𝐠_i`$ is some vector coinciding with the vector $`𝐟:=(𝚍eg(1)\mathrm{},𝚍eg(n))`$ in all the coördinates except perhaps $`i`$, and $`A`$ is the difference between the diagonal matrix $`𝚍iag\{𝐟\}`$ and the adjacency matrix of $`G`$. Clearly, $`A`$ satisfies (i)–(iv), whence $`𝐠_i`$ must  be the vector $`𝐟_i`$ obtained from $`𝐟`$ as in Section 1. Hence (3), taking in this context the form of the Coppersmith-Tetali-Winkler Identity  $`H(i,j)+H(j,k)+H(k,i)=H(j,i)+H(k,j)+H(i,k)`$, holds. Also, substituting $`𝐱_i`$ into equation $`i`$ of (6) gives the identity

$$\underset{\{k:ikG\}}{}H(k,i)=2m𝚍eg(i)$$

where $`2m:=𝚍eg(1)+\mathrm{}+𝚍eg(n)`$ is twice the number of edges in $`G`$. It follows that the return time $`R(i)`$, i.e. the expectation of the smallest natural  $`t`$ such that $`i_t=i`$ in a random walk out of $`i`$, is

$$1+\frac{1}{𝚍eg(i)}(2m𝚍eg(i))=\frac{2m}{𝚍eg(i)}.$$

Acknowledgements

I thank László Lovász, Tamás Fleiner, and Peter Winkler for their help.