The Initial Mass Function of the Inner Galaxy Measured From OGLE-III Microlensing Timescales

TL;DR

This study uses microlensing timescale data and dynamical models to measure the initial mass function in the inner Milky Way, finding it similar to local IMFs and constraining models of IMF variation.

Contribution

It provides the first measurement of the IMF in the inner Galaxy using microlensing data combined with dynamical modeling, revealing no significant difference from local IMFs.

Findings

IMF slopes consistent with Kroupa and Chabrier IMFs

Inner Galaxy IMF similar to local IMF despite old stellar population

Constraints on models of IMF variation based on gas cloud properties

Abstract

We use the timescale distribution of ~3000 microlensing events measured by the OGLE-III survey, together with accurate new made-to-measure dynamical models of the Galactic bulge/bar region, to measure the IMF in the inner Milky Way. The timescale of each event depends on the mass of the lensing object, together with the relative distances and velocities of the lens and source. The dynamical model provides statistically these distances and velocities allowing us to constrain the lens mass function, and from this to infer the IMF. Parameterising the IMF as a broken power-law, we find slopes in the main sequence $\alpha_{\rm ms}=1.31\pm0.10|_{\rm stat}\pm0.10|_{\rm sys}$ and brown dwarf region $\alpha_{\rm bd}=-0.7\pm0.9|_{\rm stat}\pm0.8|_{\rm sys}$ where we use a fiducial 50% binary fraction, and the systematic uncertainty covers the range of binary fractions 0-100%. Similarly for a log-normal IMF we conclude $M_c=(0.17\pm0.02|_{\rm stat}\pm0.01|_{\rm sys})M_\odot$ and $\sigma_m=0.49\pm0.07|_{\rm stat}\pm0.06|_{\rm sys}$. These values are very similar to a Kroupa or Chabrier IMF respectively, showing that the IMF in the bulge is indistinguishable from that measured locally, despite the lenses lying in the inner Milky Way where the stars are mostly ~10Gyr old and formed on a fast $\alpha$-element enhanced timescale. This therefore constrains models of IMF variation that depend on the properties of the collapsing gas cloud.

Figures

Tables

Table 1. Table 1 Parameters and their symmetric 1 ​ σ 1 𝜎 1\sigma errors for the power-law ( Equation 6 ) and log-normal IMF ( Equation 7 ) as a function of binary fraction. Δ ​ ln ⁡ ℒ Δ ℒ \Delta\ln\mathcal{L} is the maximum log-likelihood compared to the model with the highest log-likelihood: the log-normal IMF with 100% binary fraction.

Binary Fraction	${\alpha }_{\mathrm{ms}}$	${\alpha }_{\mathrm{bd}}$	$\Delta \ln ℒ$
0%	$1.21\pm 0.12$	$-0.02\pm 0.91$	1.1
50%	$1.31\pm 0.10$	$-0.65\pm 0.89$	0.5
100%	$1.39\pm 0.09$	$-1.45\pm 0.76$	0.4
Binary Fraction	$M_c$	${\sigma }_m$	$\Delta \ln ℒ$
0%	$0.162\pm 0.025$	$0.54\pm 0.10$	1.2
50%	$0.166\pm 0.018$	$0.49\pm 0.07$	0.5
100%	$0.169\pm 0.015$	$0.45\pm 0.05$	0.0

Equations

$$t_E=R_E/V=\frac{1}{V}\sqrt{\frac{4G{M}_l{D}_l^2}{c^2}\left(\frac{1}{D_l}-\frac{1}{D_s}\right)}$$

$$\Gamma (\log t_E|l,b,I_s,t_0)=\int$$

$$\Phi (\log M_l)\sqrt{M_l}d\log M_l$$

$${\tau }_{ij}=\frac{4\pi G}{c^2\omega }{M}_j\left(\frac{1}{D_j}-\frac{1}{D_i}\right)$$

$$\log \mathcal{L}=\sum _{\mathrm{events}i}\log {\mathcal{L}}_i.$$

$${\mathcal{L}}_i=\frac{1}{{\sigma }_{t_E,i}\sqrt{2\pi }}\int$$

$$\times \mathcal{E}(\log t_E)d\log t_E$$

$$dN$$

$$\propto M^{-\alpha }dM\text{mbox}where$$

$$\alpha$$

$$\Phi (\log M)$$

$$\alpha$$

Full text

1. Introduction

The present day mass function (PDMF), i.e. the number of stars as a function of mass, is of importance in many areas of astronomy. For example it is the key ingredient in inferring the stellar masses of galaxies from their light. The luminosity is dominated by stars with masses close to the main sequence turnoff. The mass function is then necessary to infer the total stellar mass, including the more numerous fainter dwarfs and stellar remnants. The mass distribution of stars at birth, the initial mass function (IMF), is similarly important throughout astronomy, controlling not just the PDMF but also the return of gas to, and the enrichment of, the interstellar medium. Despite the importance of the IMF we have little understanding of how it arises from the physics of the collapsing gas clouds. A variety of methods have been used to measure the PDMF and infer the IMF. The most direct measurements are from star clusters and field stars in the solar neighbourhood. Here the counting of stars makes the measurement fairly robust outside the lowest mass brown dwarfs (see e.g. Bastian et al., 2010; Krumholz, 2014; Offner et al., 2014, for recent reviews).

There are grounds for suspecting IMF variation with redshift or formation timescale because of dependence on e.g. the temperature and density of the collapsing interstellar gas (e.g. Bastian et al., 2010). There has therefore been a great effort to extend our knowledge of the IMF outside the solar neighbourhood. Mass dependent absorption features in extremely high signal to noise spectra have been used to inter the PDMF, suggesting that it may vary in massive ellipticals (e.g. van Dokkum & Conroy, 2012; Conroy & van Dokkum, 2012). Other complimentary methods to probe the IMF in external galaxies estimate dynamical masses and break the degeneracy with dark matter either through lensing (Dutton et al., 2012) or population expectations (Cappellari et al., 2012; Thomas et al., 2011).

Bulge microlensing is a unique tool for measuring the mass function in the inner Milky Way (MW), where the majority of stars formed quickly at redshift $z>1$ (i.e. they are enhanced in $\alpha$ elements and mostly $\sim 10\,{\rm Gyr}$ old, Rich, 2013). Microlensing events occur when a background source star passes in projection within the Einstein radius of a nearer star or stellar remnant and the light from the background star is therefore amplified. The level of amplification is purely geometrical and tells us no useful Galactic information. The timescale of each event however depends on several factors: the relative proper motion of the lens and source star, their distances, and the mass of the lens. Normally there is insufficient information in each event to infer all these, and so the lens mass, and the distances and velocities of the lens and source are degenerate. Dynamical models are however able to statistically provide the expected distances and velocities of microlensing events, and therefore from the distribution of microlensing timescales the lens mass distribution can be measured.

This method has been used to infer lens mass distributions several times in the literature. Both by taking moments of the distributions (Jetzer, 1994; Grenacher et al., 1999), and by fitting the timescale distribution directly from a galactic model together with a parametric IMF (Han & Gould, 1996; Zhao et al., 1996; Bissantz et al., 2004; Calchi Novati et al., 2008; Moniez et al., 2017). This work follows the latter approach.

Very recently the two most important ingredients in determining the IMF from the microlensing timescale distribution have been greatly improved, which motivates this work. Firstly Portail et al. (2017, hereafter P17) presented the first dynamical model fitted to extensive photometric and kinematic data across the bulge, bar and inner disk of the MW. This model represents a significant improvement over previously available models. Secondly Wyrzykowski et al. (2015) presented a uniform sample of 3718 microlensing events detected in Optical Gravitational Lensing Experiment (OGLE) III data, together with the necessary efficiency as a function of timescale. This is large increase over the samples of $\lesssim 50$ events used in previous microlensing IMF measurements.

2. Dynamical Model

We use the dynamical model derived from inner Galaxy data by P17 using the made-to-measure method (Syer & Tremaine, 1996; De Lorenzi et al., 2007). In this work an initially barred $N$-body model was fitted to a wide range of data consisting of: the 3D shape of the bulge measured by Wegg & Gerhard (2013), combined near-infrared star counts from the VVV, UKIDSS and 2MASS surveys (Wegg et al., 2015), and kinematics from the BRAVA (Kunder et al., 2012) and ARGOS (Ness et al., 2013) surveys. The resultant model was shown to be consistent with OGLE-II proper motion data (Sumi, 2004), important for this work which depends on these transverse velocity distributions. The reader interested in the full details of the dynamical model is referred to P17.

To provide confidence that the model reproduces the number of microlensing events we compare it to the MOA-II optical depth (Sumi & Penny, 2016) as a function of Galactic latitude in Figure 1; the model was not fitted to this data and has a reduced $\chi^{2}$ of 0.9. We also show in Figure 1 for reference how the optical depth to microlensing changes across the bulge in this model. The variation in our results from other P17 models (the boundary models), are significantly smaller than the systematic error from binaries (Section 4.1).

3. Predicted Timescale Distribution

The timescale $t_{E}$ of a microlensing event is the time for the source and lens in projection to cross the Einstein radius, $R_{E}$, of a lens with mass $M_{l}$:

[TABLE]

where $D_{l}$ and $D_{s}$ are the distances to the lens and source, and $V$ is the transverse velocity of the lens relative to the line-of-sight from Earth towards the source star.

We define $\Gamma(\log t_{E}|l,b,I_{s},t_{0})$ to be the probability distribution of $\log t_{E}$ at given position $(l,b)$, and $I$-band source magnitude $I_{s}$. The time of the event, $t_{0}$, is used to transform to the geocentric frame in which $V$ is defined.

We assume that the lens mass distribution is constant in space, and check this assumption in Section 4.2. Because $t_{E}\propto\sqrt{M_{l}}$, the timescale distribution is then given by the convolution (Han & Gould, 1996; Wegg et al., 2016)

[TABLE]

where $\gamma(\log t_{E}|l,b,I_{s},t_{0})$ is the timescale distribution predicted by the dynamical model if all lenses were $1{M_{\odot}}$, and $\Phi(\log M)$ the PDMF.

To compute $\gamma(\log t_{E}|l,b,I_{s},t_{0})$ we follow Section 3.2 of Wegg et al. (2016). In brief, at each position $(l,b)$ we select the $10^{4}$ nearest $N$-body particles as potential lens-source pairs. For each lens-source pair $(j,i)$ we compute, using Equation 1, the timescale $t_{E,ij}$ if the lens had mass $1{M_{\odot}}$. This timescale must then be correctly weighted using the $N$-body model. To do so we use the optical depth for source particle $i$ to lensing by particle $j$:

[TABLE]

where $\omega$ is the solid angle encompassed by the selected particles. This is the instantaneous probability that $j$ microlenses $i$. $\tau_{ij}M_{i}/t_{E,ij}$111This corrects a missed $j$ index from $\tau_{ij}$ in Eq (16) of Wegg et al. (2016). is therefore proportional to the expected rate of events from pair $(i,j)$. All $(i,j)$ pairs are thus binned with weight $\tau_{ij}M_{i}/t_{E,ij}$ as a function of $\log\,t_{E,ij}$ and source distance modulus $\mu_{s,i}$. Each $\mu_{s,i}$-column of the resultant matrix corresponds to the timescale distribution for that source distance modulus. However only source magnitudes are known, and not distances. We therefore convolve this with the $I$-band luminosity function of a 10$\,{\rm Gyr}$ old population to provide the rate of events at this position and time as function of source magnitude and timescale: $\gamma(\log t_{E}|l,b,t_{0},I_{s})$. Finally because here we are only interested in the distribution of timescales we normalise so that $\int\gamma(\log t_{E}|l,b,I_{s},t_{0})\,d\log t_{E}=1$.

4. Fitting the Mass Function

We utilise the timescales of the events provided by Wyrzykowski et al. (2015). We remove highly blended events with blending proportion $f_{s}<0.2$ because, as discussed there, these display a bias towards longer timescales; this leaves $2861$ events. When fitting we maximise the likelihood

[TABLE]

The likelihood of an individual event with timescale $t_{E,i}\pm\sigma_{t_{E},i}$ at $(l_{i},b_{i})$ with magnitude $I_{s,i}$ is

[TABLE]

where $\mathcal{E}(t_{E})$ is the detection efficiency of events with timescale $t_{E}$ (also provided by Wyrzykowski et al., 2015), and $\Gamma(\log t_{E}|l_{i},b_{i},{I_{s}}_{i},t_{0,i})$ is the predicted timescale distribution computed as described in the previous section. We do not focus on the shortest timescale events here, which are possibly produced by planet size mass lenses (Sumi et al., 2011; Clanton & Gaudi, 2017). We also do not fit the longest timescale events; the point lens-point source model fitted in Wyrzykowski et al. (2015) does not take account of parallax as the earth moves around the sun. These are instead the focus of dedicated searches (e.g. Wyrzykowski et al., 2016). We therefore consider only ${2\,{\rm days}}<t_{E}<200\,{\rm days}$.

We consider IMFs of broken power-law form

[TABLE]

A Kroupa (2001) IMF corresponds to $\alpha_{\rm{ms}}=1.3$ and $\alpha_{\rm{bd}}=0.3$. In the bulge $\alpha_{\rm{ms}}=1.25\pm 0.19$ was measured from star counts with an only slightly different break ($M_{\rm break}=0.56M_{\odot}$, Calamida et al. 2015), quite close to the $\alpha_{\rm{ms}}=1.43\pm 0.13$ measured in the bulge at higher latitude by Zoccali et al. (2000).

We also consider IMFs of log-normal form

[TABLE]

The Chabrier (2005) IMF has $M_{c}=0.2{M_{\odot}}$ and $\sigma_{m}=0.55$.

We transform from IMF to PDMF using a $10\,{\rm Gyr}$ population and the remnant prescription of Maraston (1998). We have checked that our conclusions are insensitive to these choices by comparing to an exponentially declining star formation rate and the remnant prescription of Percival et al. (2008).

In both Equation 6 and Equation 7 we enforce continuity of the mass function at the breaks (but not its derivative) leaving two free parameters, either $(\alpha_{\rm{bd}},\alpha_{\rm{ms}})$ or $(\log M_{c},\sigma_{m})$. Throughout we have assumed flat priors on these quantities. In Figure 2 we show the resultant likelihood contours in red of $(\alpha_{\rm{bd}},\alpha_{\rm{ms}})$ above and $(\log M_{c},\sigma_{m})$ below. The best fitting timescale distribution is compared to the data in Figure 3.

4.1. The impact of binaries

Often when measuring the IMF corrections must be made for the effect of unresolved binaries. The typical Einstein radius of microlensing events towards the bulge is $2\,{\rm a.u.}$. Binaries with this separation are often poorly fit by the single lens model, frequently presenting caustic crossings. These have been removed from the sample of Wyrzykowski et al. (2015) where only events which were well fit by the single lens model are included. Binaries with separation significantly wider than $R_{E}$ will be well resolved, however binaries with smaller separation will be unresolved. Here we assess the impact of this on the IMF. To do so we perform a binary population synthesis using the code of Hurley et al. (2002).

We evolve a population of solar metallicity binaries characterised by primary mass $M_{1}$, secondary mass $M_{2}$, separation $a$, eccentricity $e$ for 10$\,{\rm Gyr}$. From solar neighbourhood data (Raghavan et al., 2010) we assume:

(a) The period distribution $P(\log P)\propto\exp\left[\left(\log P-\log P_{0}\right)^{2}/2\sigma_{P}^{2}\right]$ with $\log(P_{0}/{\rm days})=5.0$ and $\sigma_{P}=2.3$.

(b) A flat mass ratio, $P(q\equiv M_{1}/M_{2})=\mbox{const.}$, between 0 and 1.

(c) An initially thermal eccentricity distribution, $P(e)\propto e$ between 0 and 1.

(d) The primary mass distribution $f(M_{1})$ so that the distribution of individual stellar masses matches the single stellar case. Inclusion of the distribution of secondaries with their mass distribution means this differs from choosing $f(M_{1})$ to be Equation 6 or Equation 7.

(e) Finally because close brown dwarf companions of solar type stars are extremely rare (the brown dwarf desert, Grether & Lineweaver 2006) we remove secondaries of mass $M<M_{\rm BD}=0.08{M_{\odot}}$.

We then treat binaries separated by $<3.7\sqrt{(M_{1}+M_{2})/{M_{\odot}}}\,{\rm a.u.}$ after 10$\,{\rm Gyr}$ of evolution as unresolved (where the lens is the total system) and wider binaries as resolved (where the lens is one component of the system). Since we are concerned with estimating the size of the correction more sophisticated simulations are unwarranted.

We show in Figure 2 the effect of changing the IMF to binary fractions of 50% and 100% (50% corresponding to two-thirds of all stars born in binaries). We use the locally measured 50% as our fiducial binary fraction (Raghavan et al., 2010), but because the bulk of the lens population are M-dwarfs and the binary fraction of such low-mass stars in the bulge has not been measured, we also consider 0% and 100% binary fractions.

4.2. Summary of Results

Computing confidence intervals from the likelihood assuming flat priors on the IMF parameters gives the values in Table 1. Different binary fractions and IMF forms lie within $\Delta\ln\mathcal{L}<1.2$ of each other and therefore given the uncertainties in binary population we do not significantly prefer any.

Figure 2 shows that formally the IMF differs statistically at the $\sim 2\sigma$ level with the local disk values from Kroupa (2001) and Chabrier (2005). This suggests that the average mass is lower than those fiducial values. However the differences of $\Delta\alpha_{\rm{ms}}\sim 0.1$ and $\Delta M_{c}\sim 0.02{M_{\odot}}$ are within the error budget of local IMF determinations. We therefore conclude that the inner MW IMF measured here is indistinguishable from that measured in the local disk.

Taking the locally measured 50% binary fraction as the fiducial value we conclude for the power-law IMF $\alpha_{\rm{ms}}=1.31\pm 0.10|_{\rm stat}\pm 0.10|_{\rm sys}$ and $\alpha_{\rm{bd}}=-0.7\pm 0.9|_{\rm stat}\pm 0.8|_{\rm sys}$ where the systematic uncertainty covers the range of binary fractions $0-100\%$. Similarly for the log-normal IMF we conclude $M_{c}=(0.17\pm 0.02|_{\rm stat}\pm 0.01|_{\rm sys}){M_{\odot}}$ and $\sigma_{m}=0.49\pm 0.07|_{\rm stat}\pm 0.06|_{\rm sys}$.

Note that from Figure 2 $\alpha_{\rm{bd}}<1$ at $1\sigma$ for all binary fractions meaning that the number (and the mass) per logarithmic mass interval is falling towards lower mass.

As a consistency check Figure 4 shows the variation of $\left<\log t_{E}\right>$ in Galactic coordinate bins. The trend seen in the data of shorter timescales close to the Galactic centre is reproduced by the model. Significant changes in the PDMF over this area, or differences between dynamical model and the MW, would show as inconsistent model variation.

4.3. Stellar Mass Black Holes

We show in the lower panel of Figure 3 the timescale distribution in our fiducial model, and of each type of lens separately. While the overall fit is generally excellent, for the longest events with $t_{E}\gtrsim 100\,\mbox{days}$, the fiducial model predicts more events than observed. There are several possible explanations:

Not including parallax motion when fitting the light curve in the simple lens model fitted by Wyrzykowski et al. (2015) could result in longer timescale events being inefficiently detected or $t_{E}$ misestimation. 2. 2.

This discrepancy occurs where black hole lenses become important and so could indicate the remnant prescription produces too much mass in BHs. 3. 3.

Part of the discrepancy could be because BHs and NSs are likely to receive significant natal kicks (e.g. Hansen & Phinney, 1997; Repetto et al., 2012). The distribution of remnants will therefore expand in space, particularly in the disk, away from the low-latitude sight-lines towards the bulge. In addition the remnants that do microlens would have shorter timescales because of their larger velocities. We show the effect of a Maxwellian natal kick distribution of $\sigma_{\rm{1d}}=190\,{\rm km}\,{\rm s}^{-1}$ as the dotted lines in the lower panel of Figure 3. The distribution and kinematics of NS and BHs were computed by applying an impulsive kick from this distribution to the stars in the model, and integrating forwards 1.8$\,{\rm Gyr}$. All other details of the best fitting model were kept fixed.

The other dynamical models tested from P17 give very similar long timescale behaviour (as asymptotically expected, Mao & Paczynski, 1996) and so the discrepancy is unlikely to result from the choice of model. We have also checked that reducing the upper timescale limit to $100\,{\rm days}$ does not change our conclusions, altering the best fitting model by well within the statistical errors.

5. Discussion

5.1. Comparison to other methods of IMF determination

Here we measure the IMF towards the inner MW using microlensing. The most directly comparable method is that of star counts towards the bulge. Microlensing has the advantage that it readily reaches lower masses than accessible to star counts, even down to brown dwarf and planet masses.

Local measurements using star counts in the field and nearby star clusters have uncertainties due to unresolved stellar multiplicity and, particularly for brown dwarfs, masses that can depend significantly on uncertain evolutionary models.

In comparison microlensing timescales are fundamentally sensitive directly to the mass of the lens. Conversion of timescale to mass does require a Galactic model, however because of advances in modelling techniques and the wealth of kinematic data on the inner Galaxy this uncertainty is now small. In addition only binaries with separation less than a few $\,{\rm a.u.}$ are unresolved, making inferences on the mass of the individual stars more straightforward compared to photometric IMF determinations.

5.2. The Mass-to-Light and Mass-to-Clump Ratios of the Bulge

Even in deep star counts, only stars with $M\gtrsim 0.15{M_{\odot}}$ are presently observable in the Bulge and therefore when calculating mass-to-light assumptions must be made about the contribution of lower mass stars and remnants. Microlensing is sensitive to these unseen objects and so the measurements here can be combined with star counts to robustly measure the mass-to-light. Calamida et al. (2015) measured the mass of visible stars with $>0.16{M_{\odot}}$ in the SWEEPS field to be $137,500\pm 23,400{M_{\odot}}$. These stars produced $L_{F814W}=104,000\pm 2,000L_{\odot}$ and $L_{F606W}=71,000\pm 1,400L_{\odot}$. Using the PDMF from the power-law IMFs consistent with the microlensing data gives $(M/L)_{F814W}=(2.2\pm 0.3){M_{\odot}}/L_{\odot}$ and $(M/L)_{F606W}=(3.2\pm 0.6){M_{\odot}}/L_{\odot}$ where the uncertainties are dominated by the uncertainty in the mass of the visible stars.

The dynamical models of the bulge in P17 used red clump stars (RCGs) as tracers of the MW’s structure and therefore used a mass-to-clump ratio analogous to mass-to-light in external galaxies. P17 used $(1000\pm 100){M_{\odot}}/{\rm{RCG}}$ but this required assumptions about the mass of low mass stars and remnants. It was based on $2255{M_{\odot}}/\,{\rm arcmin}^{2}$ for stars $>0.15{M_{\odot}}$ at $(l,b)=(0\deg,-6\deg)$ and $4.0\pm 0.4\mbox{RCGs}/\,{\rm arcmin}^{2}$ in the same direction (Zoccali et al., 2000). Computing the mass-to-clump similarly to the mass-to-light above gives $(960\pm 100){M_{\odot}}/{\rm{RCG}}$. Because this agrees with the mass-to-clump used by P17 it strengthens the argument made there for a low dark matter fraction in the Bulge (Wegg et al. 2016; P17) which to be reconciled with the circular velocity and dark matter estimates locally requires a core in the MW’s dark matter halo .

5.3. The Inner MW IMF in Context

The P17 model predicts that the lenses have mean distance 6.3$\,{\rm kpc}$, with a mean galactocentric radius of 2.0$\,{\rm kpc}$. As such the PDMF measured from the microlensing timescale distribution probes the IMF in the bulge and disk of the inner Galaxy. It does so directly through the lens mass, uniquely down to low masses. The much larger sample of microlensing events from OGLE-III, combined here with the newly constructed dynamical models of the inner Galaxy, provide more stringent constraints on the IMF of the inner Galaxy than previously possible.

We find that the IMF in the inner Galaxy is consistent with those measured locally by Kroupa (2001) and Chabrier (2005). However the inner Galaxy formed on a much shorter timescale than the local disk: it is $\alpha$-enhanced with a formation timescale $\sim 0.5\,{\rm Gyr}$ (Matteucci, 2014). It is also significantly older: most bulge stars are $\sim 10\,{\rm Gyr}$ (Clarkson et al., 2011; Bensby et al., 2017). The consistency of the IMF between the inner MW measured here and the local disk therefore places stringent constraints on star formation models where the IMF varies according to the properties of the parent molecular gas cloud (see e.g. Guszejnov et al., 2017, and references therein).

The authors thank the anonymous referee for their timely and thorough reading of the manuscript.