A MEASUREMENT OF DIFFUSION IN 47 TUCANAE

We inserted $\sim {{10}^{8}}$ artificial stars into the WFC3 images in F225W and F336W over the full range of observed magnitudes in both bands and a range of distances from the center of the cluster. To determine the completeness rate for the WDs that we have observed, we inserted artificial stars whose F225W and F336W magnitudes lie along the observed WD track in the CMD. The rate of recovering a star along the WD track of a given input magnitude in F336W at a given radius is the completeness rate and is depicted in Figure 2. If an artificial star along the WD track is detected in F336W, it is always detected in F225W as well. The completeness rate is both a strong function of radius and magnitude and is significantly different from unity except for the brightest stars, so accounting for completeness robustly is crucial in the subsequent analysis. The radial bins are 100 pixels in width and the magnitude bins are 0.1358 wide.

Zoom In Zoom Out Reset image size

The magnitudes of the recovered stars give the error distribution as a function of the input magnitude and position of the star in the field. Furthermore, these distributions are not typically normal and often asymmetric as well. For the analysis in Section 3.5 we use the cumulative distribution of magnitude errors as a function of position and input magnitude, which we obtain by sorting the output magnitudes in a given bin and spline to obtain the cumulative distribution in the form of the values of the errors from the first to the ninety-ninth percentile. In the analysis the completeness rate is interpolated over the two dimensions of radius and magnitude with a third-degree spline, and the error distributions are interpolated linearly over the three dimensions (radius, magnitude, and percentile).

Sometime during the late evolution of a turn-off star in 47 Tuc, the star loses about 40% of its mass, going from a main-sequence star of 90% of a solar mass to a WD of 53% of a solar mass (Renzini & Fusi Pecci 1988; Renzini et al. 1996; Moehler et al. 2004; Kalirai et al. 2009). These newborn WDs will have the typical velocities of their more massive progenitors, so as they interact gravitationally with other stars, their velocities will increase through two-body relaxation, bringing their kinetic energies into equipartition (e.g., Spitzer 1987). Because the gravitational interaction is long-range and the distance between the stars is small compared to the size of the cluster, the change in velocity will be dominated by distant interactions and small random velocity jumps, i.e., the Coulomb logarithm is large, $\sim {\rm ln} N$, where N is the number of stars in the core of the cluster, $\sim {{10}^{5-6}}$. These small jumps in velocity can be modeled as a random walk in velocity so the square of the velocity increases linearly in time and the relaxation time can be defined as ${{t}_{r}}={{\left[ d{\rm ln} {{v}^{2}}/dt \right]}^{-1}}.$ In the center of the cluster, the density of stars is approximately constant, so the gravitational potential has the approximate form

$$\begin{eqnarray}&&\phi =\frac{2}{3}\pi G{{\rho }_{c}}{{r}^{2}}.\end{eqnarray} \tag{ 1 }$$

By the virial theorem the mean kinetic energy of the WDs will equal the mean potential energy,

$$\begin{eqnarray}&&\frac{1}{2}\langle {{v}^{2}}\rangle =\frac{2}{3}\pi G{{\rho }_{c}}\langle {{r}^{2}}\rangle \end{eqnarray} \tag{ 2 }$$

so the square of the distance of the WDs from the center of the cluster will also increase linearly with time as a random walk; therefore, let us suppose that newly born WDs diffuse outward through the cluster following the diffusion equation

$$\begin{eqnarray}&&\frac{\partial \rho (r,t)}{\partial t}=\kappa {{\nabla }^{2}}\rho (r,t),\end{eqnarray} \tag{ 3 }$$

where κ can be related to the relaxation time as $\kappa =r_{c}^{2}/{{t}_{r}}$ because $d{\rm ln} {{r}^{2}}/dt=d{\rm ln} {{v}^{2}}/dt$ from Equation (2). This diffusion equation yields the Green's function

$$\begin{eqnarray}&&u(r,t)=\frac{1}{8{{\left( \pi \kappa t \right)}^{3/2}}}{{e}^{-{{r}^{2}}/(4\kappa t)}}\end{eqnarray} \tag{ 4 }$$

if κ is independent of time and position. This gives a cumulative distribution in projected radius

$$\begin{eqnarray}&&C(\lt R)=1-{{e}^{-{{R}^{2}}/(4\kappa t)}}.\end{eqnarray} \tag{ 5 }$$

The Green's function at t = 0 is a delta function centered on the center of the cluster. On the other hand, if the initial distribution is a Gaussian centered on the center of the cluster the density of WD stars near the center of the cluster is a function of age, t, and projected radius, R, of the form

$$\begin{eqnarray}&&\rho =\rho (R,t)=\frac{2R}{4\kappa (t+{{t}_{0}})}{\rm exp} \left[ -\frac{{{R}^{2}}}{4\kappa (t+{{t}_{0}})} \right]\end{eqnarray} \tag{ 6 }$$

where the density distribution is normalized as

$$\begin{eqnarray}&&\int _{0}^{\infty }\rho (R,t)dR=1.\end{eqnarray} \tag{ 7 }$$

The dispersion of the Gaussian at t = 0 is simply given by $2{{\sigma }^{2}}=4\kappa {{t}_{0}}$. Because the diffusion equation is linear, a sum of several Gaussians with the same value of κ but different normalizations and values of t₀ will also be a solution.

Of course we do not directly observe the ages of the WDs. Rather we observe their fluxes or apparent magnitudes. The cooling curve of the WDs is a relationship between time and the apparent magnitude from the WDs t(m), so the number of WDs that we expect to observe at a given flux and radius is given by

$$\begin{eqnarray}&&f(R,m)=\dot{N}\rho (R,t(m))\frac{\partial t}{\partial m}C(R,m),\end{eqnarray} \tag{ 8 }$$

where $\dot{N}$ is the birthrate of the WDs (assumed to be constant over the range of ages of the young WDs, i.e., the past 200 Myr) and $C(R,m)$ is the completeness as a function of radius and flux. To this point flux errors have been neglected.

To construct the various cooling models here, i.e., t(m) from Equation (8), we used Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011) to perform simulations of stellar evolution starting with a pre-main-sequence model of 0.9 solar masses and a metallicity of $Z=3.3\times {{10}^{-3}}$ appropriate for the cluster 47 Tuc. This is slightly larger than the value for the turnoff mass found by Thompson et al. (2010) for the eclipsing binary V69 in 47 Tuc that is composed of an upper main sequence star of about 0.86 solar masses and a subgiant of 0.88 solar masses. Because we are interested in the stars that have become young WDs just recently, the initial masses of these stars should be slightly larger than the turnoff mass today. We have explicitly assumed that the progenitors of the WDs are a uniform population. Although there is evidence of modest variation in the chemical abundances in 47 Tuc (e.g., Milone et al. 2012), the WD cooling sequence, at least at larger radii, appears uniform (Richer et al. 2013). However, from Figure 1 it is apparent that the core of 47 Tuc has a substantial population of blue stragglers that will evolve to become more massive WDs. Our sample has about 160 blue stragglers, and if we estimate the duration of the main-sequence for a blue straggler to be 1–2 Gyr (e.g., Sills et al. 2009), we obtain a birth-rate of blue-straggler WDs of about 0.1 Myr⁻¹. The number of giants in our field indicates a birth rate of about eight WDs per million years (see Goldsbury et al. 2012 for further details), so the estimated contamination of the WD cooling sequence is modest at about 1%.

Specifically, we used SVN revision 5456 of MESA and started with the model 1 M_pre_ms_to_wd in the test suite. We changed the parameters initial_mass and initial_z of the star and adjusted the parameter log_L_lower_limit to −6 so the simulation would run well into the WD cooling regime. We also reduced the two values of the wind η to 0.46 (from the default of 0.7) to yield a 0.53 solar mass WD from the 0.9 solar mass progenitor. Interestingly Miglio et al. (2012) argue from Kepler asteroseismic measurements of the stars in the metal-rich open cluster NGC 6971 that such values of η are needed to account for the mass loss between the red giant and red clump phases of stars in this metal-rich cluster.

We defined the time of birth of the WD to coincide with the peak luminosity of the model at the tip of the asymptotic giant branch about 10.9 Gyr after the start of the simulation. This is in agreement with the best age of the cluster determined from main-sequence stars of 11.25 ± 0.21 (random) ±0.85 (systematic) Gyr (Thompson et al. 2010). This age agrees with that derived by Hansen et al. (2013) from WD cooling (9.9 ± 0.7 Gyr at 95% confidence). We choose this definition of the birth so that each observed WD will have a star of similar luminosity in the cooling model. At this point in the evolution we have outputs from the MESA evolution every 100 yr or so; therefore, the cooling curve is well-sampled throughout. At each output time we have the value of the luminosity, radius, effective temperature, and mass of the star. With these values we interpolate the spectral models of Tremblay et al. (2011) in surface gravity and effective temperature, and then scale the result to the radius of the model star. We use a true distance modulus of 13.23 (Thompson et al. 2010) and a reddening of $E(B-V)=0.04$ (Salaris et al. 2007) to determine the model fluxes in the WFC3 band F336W. We used the standard extinction curve of Fitzpatrick (1999) with ${{A}_{V}}/E(B-V)=3.1$. We have purposefully used a distance and reddening determined from main-sequence stars to avoid a potential circularity in using the WD models themselves to fix the distance. Woodley et al. (2012) inferred a slightly larger true distance modulus of $13.36\;\pm \;0.02\;\pm \;0.06$ from the WD spectral energy distributions.

The brightest WD in our sample has F336W = 14.92. According to the models this corresponds to an age of 110,000 yr, an effective temperature of 100,000 K, a luminosity of 1600 L_⊙, and a radius of 0.13 R_⊙. Its mass is 0.53 solar masses. The faintest WD in our sample has F336W = 25.4, yielding an age of 1.2 Gyr, an effective temperature of 8700 K, a luminosity of 10⁻³ L_⊙, and a radius of 0.013 R_⊙, one-tenth of the radius of the brightest WD. Clearly the brightest WD in our sample is not a WD in the usual sense because thermal energy plays an important role in the pressure balance of the star. For this brightest star ${\rm log} g=5.93$, which is less than the minimum of the atmosphere model grid (${\rm log} g=6$) so we have to extrapolate slightly off of the grid, but only for this brightest star. For the simulations in Section 3.8 we did not use this particular model, but similar ones of the same WD mass with different neutrino cooling rates or initial metallicities also generated with MESA.

The model outlined in Section 3.2 predicts the number of WDs as a function of magnitude and position. Let us divide the space of position and magnitude into bins of width ${\Delta }R$ and ${\Delta }m$ where the bins are numbered with indices j and k, respectively. The probability of finding n stars in a particular bin is given by

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} P(n;f({{R}_{j}},{{m}_{k}})) \\ \quad =\frac{{{\left[ f({{R}_{j}},{{m}_{k}}){\Delta }R{\Delta }m \right]}^{n}}{{e}^{-f({{R}_{j}},{{m}_{k}}){\Delta }R{\Delta }m}}}{n!}. \\ \end{array}\end{eqnarray} \tag{ 9 }$$

Now we imagine dividing the sample into so many bins that there is either a single star in a bin or no stars at all; we have

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} P(n;f({{R}_{j}},{{m}_{k}})) & = & {{e}^{-f({{R}_{j}},{{m}_{k}}){\Delta }R{\Delta }m}} \\ {} & {} & \times \left\{ \begin{array}{ccccccccccccccc} f({{R}_{j}},{{m}_{k}}){\Delta }R{\Delta }m & {\rm one}\;{\rm star} \\ 1 & {\rm no}\;{\rm star} \\ \end{array} \right.. \\ \end{array}\end{eqnarray} \tag{ 10 }$$

We can define the likelihood as the logarithm of the product of the probabilities of observing the number of stars in each bin. Since the bins are so small we can replace R_j and m_k for the bins with stars in them with the measured values for that particular star R_i and m_i. This gives the so-called "unbinned likelihood" of observing the sample as follows (Cash 1979; Mattox et al. 1996; Davis 2014):

$$\begin{eqnarray}&&{\rm log} L=\mathop{\sum }\limits_{i}{\rm log} f({{R}_{i}},{{m}_{i}})-\mathop{\sum }\limits_{j,k}f({{R}_{j}},{{m}_{k}}){\Delta }R{\Delta }m.\end{eqnarray} \tag{ 11 }$$

We have dropped the constant widths of the bins from the first term, which is a sum over the observed stars; consequently, the absolute value of the likelihood is not important, just differences matter. The second term is a sum over the really narrow (and arbitrary) bins that we have defined, so we have

$$\begin{eqnarray}&&\mathop{\sum }\limits_{j,k}\;f({{R}_{j}},{{m}_{k}}){\Delta }R{\Delta }m=\int f({{R}_{j}},{{m}_{k}})dRdm={{N}_{{\rm pred}}},\end{eqnarray} \tag{ 12 }$$

where N_pred is the number of stars that the model predicts that we will observe, so finally we have

$$\begin{eqnarray}&&{\rm log} L=\mathop{\sum }\limits_{i}{\rm log} f({{R}_{i}},{{m}_{i}})-{{N}_{{\rm pred}}},\end{eqnarray} \tag{ 13 }$$

where the summation is over the observed stars. The integral for N_pred when combined with Equation (8) yields

$$\begin{eqnarray}&&{{N}_{{\rm pred}}}=\dot{N}\int _{0}^{{{t}_{1}}}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,t)C(R,m)dRdt\end{eqnarray} \tag{ 14 }$$

or

$$\begin{eqnarray}&&{{N}_{{\rm pred}}}=\dot{N}\int _{{{m}_{0}}}^{{{m}_{1}}}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,t)C(R,m)\frac{dt}{dm}dRdm.\end{eqnarray} \tag{ 15 }$$

If we take the luminosity function as fixed and try to maximize the likelihood with respect to the diffusion model

$$\begin{eqnarray}&&{\rm log} L=\mathop{\sum }\limits_{i}{\rm log} \left[ \dot{N}\rho (R,t({{m}_{i}})){{\left. \frac{\partial t}{\partial m} \right|}_{i}}C({{R}_{i}},{{m}_{i}}) \right]-{{N}_{{\rm pred}}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {} & = & \mathop{\sum }\limits_{i}{\rm log} \left[ \dot{N}\rho (R,t({{m}_{i}})) \right] \\ {} & {} & +\mathop{\sum }\limits_{i}{\rm log} \left[ {{\left. \frac{\partial t}{\partial m} \right|}_{i}}C({{R}_{i}},{{m}_{i}}) \right]-{{N}_{{\rm pred}}}, \\ \end{array}\end{eqnarray} \tag{ 17 }$$

where the second summation does not depend on the diffusion model, so it is constant with respect to changes in the diffusion model and can be dropped from the logarithm of the likelihood. However, it must be included if one wants to compare different cooling curves, t(m).

An important complication to the analysis is that the measured magnitudes are not the same as the actual magnitudes of the stars; in particular the error distribution is not normal or even symmetric. This transforms the model distribution function via a convolution,

$$\begin{eqnarray}&&f^{\prime} (R,m)=\int _{-\infty }^{\infty }f(R,m^{\prime} )g(R,m^{\prime} ,m-m^{\prime} )dm^{\prime} \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {} & = & \int _{-\infty }^{\infty }f(R,m-{\Delta }m)g(R,m-{\Delta }m,{\Delta }m)d\left( {\Delta }m \right) \\ {} & = & \int _{0}^{1}f(R,m-{\Delta }m)dG, \\ \end{array}\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}&&G(R,m,{\Delta }m)=\int _{-\infty }^{{\Delta }m}g(R,m-{\Delta }m,{\Delta }m^{\prime} )d\left( {\Delta }m^{\prime} \right)\end{eqnarray} \tag{ 20 }$$

is the cumulative distribution of magnitude errors with the observed radius and magnitude fixed. If we calculate the percentiles of the magnitude error distribution as ${\Delta }{{m}_{j}}$ we can approximate the integral as the sum

$$\begin{eqnarray}&&f^{\prime} (R,m)=\frac{1}{100}\mathop{\sum }\limits_{j}\;f(R,m-{\Delta }{{m}_{j}}),\end{eqnarray} \tag{ 21 }$$

so for a given star i we have

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} f^{\prime} ({{R}_{i}},{{m}_{i}}) & = & \dot{N}\mathop{\sum }\limits_{j}\rho (R,t({{m}_{i}}-{\Delta }{{m}_{j}})) \\ {} & {} & \times t^{\prime} \left( {{m}_{i}}-{\Delta }{{m}_{j}} \right)C({{R}_{i}},{{m}_{i}}-{\Delta }{{m}_{j}}), \\ \end{array}\end{eqnarray} \tag{ 22 }$$

where $t^{\prime} =\partial t/\partial m$. This new function $f^{\prime} ({{R}_{i}},{{m}_{i}})$ can be substituted into Equation (13) to yield a likelihood including the magnitude errors. We will assume that the magnitude errors do not affect our estimate of N_pred; this simplifies the analysis. We will verify our technique with Monte Carlo simulations in Section 3.8.

We can construct a maximum likelihood estimator of the luminosity function of the WDs or alternatively the cooling curve as follows

$$\begin{eqnarray}&&f(R,m)=\frac{dN}{dmdR}=\dot{N}\rho (R,t)\frac{\partial t}{\partial m}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&=\;\dot{N}\rho (R,t)\mathop{\sum }\limits_{i}{{A}_{i}}\delta (m-{{m}_{i}}),\end{eqnarray} \tag{ 24 }$$

where i is an index that runs over the observed stars. With this model we can define a likelihood function for the stars that we observe

$$\begin{eqnarray}&&{\rm log} L=\mathop{\sum }\limits_{i}{\rm log} \left[ \dot{N}{{A}_{i}}\rho ({{R}_{i}},{{t}_{i}}) \right]-{{N}_{{\rm pred}}},\end{eqnarray} \tag{ 25 }$$

where a multiplicative constant (infinite in this case) and the completeness for each star have been dropped from the logarithm.

Substituting the trial luminosity function Equation (23) yields

$$\begin{eqnarray}&&{{N}_{{\rm pred}}}=\mathop{\sum }\limits_{i}{{A}_{i}}\dot{N}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,{{t}_{i}})C(R,{{m}_{i}})dR.\end{eqnarray} \tag{ 26 }$$

If we maximize the likelihood with respect to the values of A_k we obtain

$$\begin{eqnarray}&&\frac{\partial {\rm log} \left[ \rho ({{R}_{i}},{{t}_{i}})C({{R}_{i}},{{m}_{i}}) \right]}{\partial {{A}_{k}}}=\frac{\partial {\rm log} \left[ \rho ({{R}_{i}},{{t}_{i}}) \right]}{\partial {{t}_{i}}}\frac{\partial {{t}_{i}}}{\partial {{A}_{k}}},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}\frac{\partial {{t}_{i}}}{\partial {{A}_{k}}}=\left\{ \begin{array}{ccccccccccccccc} 1 & {\rm if} & k\leqslant i \\ 0 & {\rm if} & k\gt i \\ \end{array} \right.\end{eqnarray} \tag{ 28 }$$

so

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} \mathop{\sum }\limits_{i}\frac{\partial {\rm log} \left[ {{A}_{i}}\rho ({{R}_{i}},{{t}_{i}})C({{R}_{i}},{{m}_{i}}) \right]}{\partial {{A}_{k}}} \\ \quad =\frac{1}{{{A}_{k}}}+\mathop{\sum }\limits_{i\geqslant k}\frac{\partial {\rm log} \left[ \rho ({{R}_{i}},{{t}_{i}}) \right]}{\partial {{t}_{i}}} \\ \end{array}\end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&\frac{\partial {\rm log} \left[ \rho (R,t) \right]}{\partial t}=\frac{{{R}^{2}}}{4\kappa {{(t+{{t}_{0}})}^{2}}}-\frac{1}{t+{{t}_{0}}}.\end{eqnarray} \tag{ 30 }$$

Taking the derivative of N_pred yields the second part of the variance in ${\rm log} L$,

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{\partial {{N}_{{\rm pred}}}}{\partial {{A}_{k}}} & = & \dot{N}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,{{t}_{k}})C(R,{{m}_{k}})dR \\ {} & {} & +\mathop{\sum }\limits_{i}{{A}_{i}}\dot{N}\int _{0}^{{{r}_{{\rm max} }}}\frac{\partial \rho (R,{{t}_{i}})}{\partial {{t}_{i}}}\frac{\partial {{t}_{i}}}{\partial {{A}_{k}}}C(R,{{m}_{i}})dR \\ \end{array}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {} & = & \dot{N}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,{{t}_{k}})C(R,{{m}_{k}})dR \\ {} & {} & +\mathop{\sum }\limits_{i\geqslant k}{{A}_{i}}\dot{N}\int _{0}^{{{r}_{{\rm max} }}}\frac{\partial \rho (R,{{t}_{i}})}{\partial {{t}_{i}}}C(R,{{m}_{i}})dR. \\ \end{array}\end{eqnarray} \tag{ 32 }$$

Combining these results with $\partial {\rm log} L/\partial {{A}_{k}}=0$ yields an equation of the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{1}{{{A}_{k}}} & = & \dot{N}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,{{t}_{k}})C(R,{{m}_{k}})dR \\ {} & {} & \mathop{\sum }\limits_{i\geqslant k}\left[ {{A}_{i}}\dot{N}\int _{0}^{{{r}_{{\rm max} }}}\frac{\partial \rho (R,{{t}_{i}})}{\partial {{t}_{i}}}C(R,{{m}_{i}})dR \right. \\ {} & {} & -\left. \frac{\partial {\rm log} \left[ \rho ({{R}_{i}},{{t}_{i}}) \right]}{\partial {{t}_{i}}} \right] \\ \end{array}\end{eqnarray} \tag{ 33 }$$

or a matrix equation of the form

$$\begin{eqnarray}&&\mathop{\sum }\limits_{i}{{M}_{ik}}{{A}_{i}}={{b}_{k}}+\frac{1}{{{A}_{k}}},\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}{{M}_{ik}}=\left\{ \begin{array}{ccccccccccccccc} {{M}_{i}} & {\rm if} & i\geqslant k \\ 0 & {\rm if} & i\lt k \\ \end{array} \right.\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&{{M}_{i}}=\dot{N}\int _{0}^{{{r}_{{\rm max} }}}\frac{\partial \rho (R,{{t}_{i}})}{\partial {{t}_{i}}}C(R,{{m}_{i}})dR.\end{eqnarray} \tag{ 36 }$$

The vector b_k is given by

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{b}_{k}} & = & \mathop{\sum }\limits_{i\geqslant k}\frac{\partial {\rm log} \left[ \rho ({{R}_{i}},{{t}_{i}}) \right]}{\partial {{t}_{i}}}-\dot{N} \\ {} & {} & \times \int _{0}^{{{r}_{{\rm max} }}}\rho (R,{{t}_{k}})C(R,{{m}_{k}})dR. \\ \end{array}\end{eqnarray} \tag{ 37 }$$

Although this matrix equation has as many rows as there are stars in the sample, it is straightforward to solve at least formally in two ways. The first is

$$\begin{eqnarray}&&{{A}_{k}}=\frac{1}{\mathop{\sum }\limits_{i}{{M}_{ik}}{{A}_{i}}-{{b}_{k}}}\end{eqnarray} \tag{ 38 }$$

and the second is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{A}_{i}} & = & \frac{{{b}_{i}}+{{({{A}_{i}})}^{-1}}}{{{M}_{i}}}-\frac{{{b}_{i+1}}+{{({{A}_{i+1}})}^{-1}}}{{{M}_{i+1}}}, \\ {{A}_{n}} & = & \frac{{{b}_{n}}+{{({{A}_{i}})}^{-1}}}{{{M}_{n}}}. \\ \end{array}\end{eqnarray} \tag{ 39 }$$

The values of b_i and M_i depend on the values of A_i through the parameter t_k, so the solution must proceed iteratively while minimizing with respect to the other parameters of the model κ and t₀.

For each value of the diffusion parameters, we chose to iterate Equation (38) three times to determine the values of A_k within a loop of two iterations where M_i (Equation (36)) and b_k (Equation (37)) vary. Given this new trial luminosity function, the diffusion parameters are varied to find the maximum likelihood, and the iterative solution of the luminosity function is repeated. These two steps are repeated until the values of the diffusion parameters from one iteration to the next have changed by less than one part per hundred.

An interesting limit is when the density distribution is independent of time. This understandably yields a simpler solution for A_i. In particular, M_i = 0 so

$$\begin{eqnarray}&&{{A}_{i}}=-\frac{1}{{{b}_{i}}}={{\left[ \dot{N}\int _{0}^{{{r}_{{\rm max} }}}\rho (R,{{m}_{i}})C(R,{{m}_{i}})dR \right]}^{-1}},\end{eqnarray} \tag{ 40 }$$

where the underlying density distribution is normalized. The weight is not the reciprocal of the completeness for star i but rather the reciprocal of the mean of the completeness of a star with the flux of star i over the density distribution. The latter could be evaluated by taking the mean of the completeness measured for all the stars in the sample in a magnitude range about star i sufficiently wide to sample the density distribution. It is important to note that the weight is the reciprocal of the mean of the completeness not the mean of the reciprocal. If the completeness does not depend strongly on radius, these two will approximately coincide. Finally if the density distribution is not known a priori and is not modeled, the weight for a particular star is simply given by ${{A}_{i}}={{\left[ C({{R}_{i}},{{m}_{i}}) \right]}^{-1}}$. We call this "Inv Comp" in Figures 6 and 10.

The likelihood is invariant under changes in the birth rate of the WDs ($\dot{N}$) if one also changes the values of $\kappa ,{{A}_{i}}$ and t₀ as follows:

$$\begin{eqnarray}&&\dot{N}\to \alpha \dot{N},{{A}_{i}}\to {{\alpha }^{-1}}{{A}_{i}},{{t}_{0}}\to {{\alpha }^{-1}}{{t}_{0}}\;{\rm and}\;\kappa \to \alpha \kappa .\end{eqnarray} \tag{ 41 }$$

That is, the timescale cannot be fixed without some additional input such as a theoretical cooling curve or an independent estimate of the WD birthrate. The quantities ${{A}_{i}}\dot{N}$, $\kappa /\dot{N}$, and $\dot{N}{{t}_{0}}$ are invariant with respect to this transformation. In our data set when we use this modeling technique, we fix the value of $\dot{N}$ to the value inferred by the number of giants in our field as in Goldsbury et al. (2012).

We start the analysis including magnitude errors with Equations (18) and (24), which when combined yield

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} f^{\prime} (R,m) & = & \int _{-\infty }^{\infty }\dot{N}\rho (R,t)C(R,m^{\prime} )\mathop{\sum }\limits_{i}{{A}_{i}}\delta (m^{\prime} -{{m}_{i}}) \\ {} & {} & \times g(R,m^{\prime} ,m-m^{\prime} )dm^{\prime} \\ \end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&=\;\dot{N}\mathop{\sum }\limits_{i}{{A}_{i}}\rho (R,{{t}_{i}})g(R,{{m}_{i}},m-{{m}_{i}}).\end{eqnarray} \tag{ 43 }$$

With this model we can define a likelihood function for the stars that we observe

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\rm log} L & = & \mathop{\sum }\limits_{j}{\rm log} \left[ \dot{N}\mathop{\sum }\limits_{i}{{A}_{i}}\rho ({{R}_{j}},{{t}_{i}})g({{R}_{j}},{{m}_{i}},{{m}_{j}}-{{m}_{i}}) \right. \\ {} & {} & \times \left. C({{R}_{j}},{{m}_{i}}) \right]-{{N}_{{\rm pred}}}. \\ \end{array}\end{eqnarray} \tag{ 44 }$$

Note how the magnitude error essentially translates into a spread in the age of the observed stars

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{\partial {\rm log} L}{\partial {{A}_{k}}} & = & \mathop{\sum }\limits_{j,i\geqslant k}\left\{ \frac{\rho ({{R}_{j}},{{t}_{i}})g({{R}_{j}},{{m}_{i}},{{m}_{j}}-{{m}_{i}})C({{R}_{j}},{{m}_{i}})}{{{\sum }_{l}}{{A}_{l}}\rho ({{R}_{j}},{{t}_{l}})g({{R}_{j}},{{m}_{l}},{{m}_{j}}-{{m}_{l}})C({{R}_{j}},{{m}_{l}})} \right. \\ {} & {} & \times \left. \left[ {{\delta }_{ik}}+{{A}_{i}}\frac{\partial {\rm log} \left[ \rho ({{R}_{i}},{{t}_{i}}) \right]}{\partial {{t}_{i}}} \right] \right\}-\frac{\partial {{N}_{{\rm pred}}}}{\partial {{A}_{k}}}. \\ \end{array}\end{eqnarray} \tag{ 45 }$$

Although we have included this additional complication in the derivations for completeness, we have found that the inclusion of error convolution in modeling simulated data does not affect the fitting results, so we did not include this in the modeling of the dynamics while simultaneously determining the luminosity functions.

To test these techniques in the face of the challenges of incompleteness and magnitude errors present in our data, we typically simulated on the order of 10,000 catalogs of the same size as our data set with a known luminosity function and a known diffusion model and attempted to recover the input parameters. In both cases, the age of the star is first selected to be between zero and 1.5 Gyr. Given this age the model cooling curve determines the F336W magnitude. Second, a radius is selected from the cumulative distribution in projected radius (Equation (5)). Given the radius and magnitude of the candidate for the catalog, the completeness for this star is calculated and the star is included in the sample with this probability. Finally, the magnitude errors are applied by drawing from the magnitude error distribution. We created a sample of 3167 stars—the same as in the WFC3 WD sample. The fitting procedure followed two different strategies.

The first was to assume a fixed cooling curve and try to find the density evolution to determine whether the process is biased in determining the diffusion parameters and the typical errors. Finally, we performed simulations where we did not convolve the models with the error distribution to calculate the likelihood (in all cases errors were applied to the simulated data) to see whether the omission of this step introduced biases. The second strategy did not assume a cooling curve and determined the cooling curve as a part of the process of determining the diffusion. We did not convolve the cooling curve with the error distribution while fitting the model; however, the fake catalogs were created in the same way as in the first strategy. In this technique the resulting cooling curve can be multiplied by a constant factor (Equation (41)), so we determine the values of κ and t₀ by fixing the value of $\dot{N}$ to the one used to build the catalog. This also fixes the age estimates of all of the WDs in the sample.

The results of these simulations are depicted in Figures 3 and 4 and in Table 1. The key results of the simulations are that the likelihood fitting of the diffusion model results in an unbiased estimate of the diffusion parameters regardless of whether the fitting technique includes the magnitude errors (Section 3.5). Furthermore, even when one fits for the luminosity function as well one can obtain reliable estimates of the diffusion model without prior knowledge of the cooling curve; of course, in this latter case the timescales of the diffusion rely on an independent estimate of the birth rate of the WDs $\dot{N}$. Observationally, this is determined from a sample of giant stars numbering in the thousands (see Figure 1) so the statistical error in this determination is small. Typically the birth rate is recovered with an uncertainty of less than 1% and the diffusion rate with an uncertainty of 10% and t₀ with an uncertainty of 15%. The errors in t₀ and κ are correlated so the error in ${{t}_{0}}\kappa$ is typically less than 10%.

Zoom In Zoom Out Reset image size

Table 1.  Parameters from the Monte Carlo Simulations: Means, Standard Derivations, Input Values

Technique	κ	Input	t0	Input
Full Modeling	3.4 ± 0.3	3.58	560 ± 80	531
No-error Convolution	3.7 ± 0.4	3.71	530 ± 80	515
Unfixed LF	7.5 ± 0.6	7.26	220 ± 30	231
	$\dot{N}$	Input	${\rm log} L$
Full Modeling	5.42 ± 0.03	5.44	±50
No-error Convolution	5.45 ± 0.03	5.45	±22
Unfixed LF	⋯	±20	⋯

Download table as:  ASCIITypeset image

In the second type of simulation, we found the density evolution along with an estimate of the cooling curve, so this cooling curve can be compared with the input cooling curve for the simulations. Furthermore, the determination of the cooling curve is iterative, so we have to give an initial guess of the curve. The input, the initial guess and the results are given in Figure 5. We can also fit for just the cooling curve and assume that the density distribution does not evolve or not assume a density model at all and use the per star completeness as outlined in Section 3.6. Figure 6 highlights the difference between the model age and the inferred age with the various likelihood techniques. For young WDs the uncertainties are large (because there are few young WDs in the sample), but for old WDs there is a small bias of order of 10% in the inferred age, the sign of which depends on the technique. Again this is on the order of the relative errors in the diffusion parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 depicts the radial distribution of WDs of various ages. Each bin is 50 Myr wide, and the bins are centered on 25, 125, and 225 Myr. The evolution at up to a few core radii (about 60 arcsec) is dramatic from 25 to 125 Myr and modest thereafter. Outside 60 arcsec the cumulative distributions are nearly parallel indicating little evolution in this region at early times. The simple diffusion models used here assume that the diffusion coefficient is constant in space and in time, so the models continue to evolve at late times and for all radii. At the smaller radii the WDs reach the distribution corresponding to their masses after about 100 Myr and stop diffusing.

Zoom In Zoom Out Reset image size

Figure 9 focuses on the outer half of the WFC3 field. Here we see more evolution between the second and third epochs, with little early evolution. This indicates the increase in the relaxation time as the stellar density decreases. The WDs diffuse modestly over the first 100 Myr and more dramatically during the second 200 Myr. The WDs as expected from theoretical considerations suffer diffusion that is a function of radius and time and beyond the scope of the simple model used to quantify the diffusion in this paper. However, this model does capture the diffusion within a few core radii for a few core relaxation times.

Zoom In Zoom Out Reset image size

In this paper we used the Green's function (Equation (4)) to model the diffusion of the stars through the cluster. We simply took the initial conditions to be a Gaussian or a sum of Gaussians centered on the center of the cluster. This allowed for a simple closed-form expression for the density function in spherical coordinates and in projection as well. Without relaxing the spherical symmetry one could imagine much more general initial conditions. In fact we have an estimate of the initial conditions in the form of the projected radial distribution of the stars on the upper main-sequence. This distribution could be possibly deprojected as a lowered-isothermal distribution in phase space (Michie 1963; King 1966) and convolved with the Gaussian Green's function, Equation (4), to give the expected density distribution as a function of time. This technique shares the advantage of the technique used in this paper that the density distribution can be guaranteed to be positive because the convolution of the positive kernel with a positive distribution is necessarily positive; however, the density distribution even in spherical coordinates is not available in closed form.

A second strategy would be to expand the initial density distribution in terms of spherical Bessel functions and spherical harmonics. If we restrict ourselves to an initially spherical distribution we have

$$\begin{eqnarray}&&\rho (r,t)=\int _{0}^{\infty }dka(k){{e}^{-\kappa {{k}^{2}}t}}\frac{{\rm sin} (kr)}{kr},\end{eqnarray} \tag{ 46 }$$

where the coefficients a(k) are determined from the initial density distribution

$$\begin{eqnarray}&&a(k)=\frac{2}{\pi }\int _{0}^{\infty }dr{{\left( kr \right)}^{2}}\rho (r,0)\frac{{\rm sin} (kr)}{kr}.\end{eqnarray} \tag{ 47 }$$

If the initial density distribution can be well represented with a few values of k, then the density evolution is straightforward to evolve forward and backward in time; however, it is no longer guaranteed to be positive even at the initial time if only a range of values of k are considered in a (k).

From the point of view of the likelihood analysis, a natural next step would be to use the additional information available with the current observations, i.e., the flux in the F225W band. This would provide an additional constraint on the ages of the WD stars or alternatively constrain the cooling curve in both bands. In the first case one would perhaps get better constraints on the dynamical evolution and could also fit for the distance and reddening to the cluster, and possibly the mass of the WDs or specifics of the cooling mechanism. In the second case one would get a cooling curve in a second band. It is straightforward to see that the weights for the cooling curve in F225W would be the same as in F336W, so simply plotting the inferred ages of the WDs from Figure 10 against the F225W magnitude would yield the cooling curve in F225W. The agreement with the F336W model is poorer at early times, but improves with age and lasts until nearly 1 Gyr. In the context of this paper, we obtain similar diffusion parameters whether we fit a luminosity function or assume a theoretical model.

Zoom In Zoom Out Reset image size

As argued in Section 3.2 the interactions with other stars cause the WDs to diffuse in velocity, not in radius. However, using the virial theorem we argued that this diffusion in velocity would be manifest as a diffusion in radius as well. Furthermore, our simple model assumes that the diffusion coefficient is constant in space and time when in fact with time the WD distribution approaches that of stars of similar mass so the diffusion must cease and also at larger radii the diffusion must happen more slowly. We see both of these effects in Figures 8 and 9. How this diffusion actually manifests itself could be simulated in two possible ways.

The first is direct numerical simulation of on order of one million stars that form the central regions of the globular cluster 47 Tuc. This appears to be a Herculean labor because the state-of-the-art direct calculation of the two body interactions in a globular cluster involve merely $\sim {{10}^{5}}$ stars and the simulation in question would normally take 100 times longer. However, we are only interested in the dynamical evolution of the young WDs over about one hundredth of the age of the cluster (100–150 Myr out of 10 Gyr). Second, because we are not interested in the long-term evolution of the cluster, neither stellar evolution nor the dynamics of binaries should play an important role in this process. These two simplifications result in a factor of a thousand speed up to obtain results and these calculations are already underway.

The second direction would be to model the diffusion in phase space using a Fokker–Planck or Monte Carlo scheme (Giersz & Heggie 2011; Pattabiraman et al. 2013; Hong et al. 2013). Such simulations would be more rapid than a direct N-body simulations and possibly yield more physical insights.

Following the arguments of Section 5.2 a natural direction would be to measure the proper motions of the WDs in the core of 47 Tuc with a second epoch of observations. Because we already have the colors of the WDs, only observations in a single band would be required and possibly not as deep as the present set of observations because the stars have already been detected. To obtain the most precise positions and to minimize the crowding, the bluest band would be best, i.e., F225W, and possibly over only a portion of the field of the current data, because here the goal would be to verify the current result by finding the corresponding signal in velocity space, so a full sample of 3000 plus WDs may not be required.

We have measured directly for the first time the dynamical relaxation of stars in a globular cluster. To do this we have introduced new statistical techniques for the characterization of stellar populations. These techniques can robustly and straightforwardly account for high incompleteness and non-Gaussian magnitude errors. They can be applied to a wide variety of questions from globular cluster dynamics to galaxy luminosity functions. There are many avenues for further investigation, such as a more thorough analysis of the existing data using the information from the second band, the simulation of the relaxation of young WDs in numerical models, and measuring the proper motions of the young WDs to search for signatures of relaxation in their velocities as well.

This research is based on NASA/ESA Hubble Space Telescope observations obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy Inc. under NASA contract NAS5-26555. These observations are associated with proposal GO-12971 (PI: Richer). This work was supported by NASA/HST grant GO-12971, the Natural Sciences and Engineering Research Council of Canada, the Canadian Foundation for Innovation, and the British Columbia Knowledge Development Fund. This project was supported by the National Science Foundation (NSF) through grant AST-1211719. It has made used of the NASA ADS and arXiv.org.