RECONSTRUCTING THE ACCRETION HISTORY OF THE GALACTIC STELLAR HALO FROM CHEMICAL ABUNDANCE RATIO DISTRIBUTIONS

The simulations consist of 11 "MW-sized" halo realizations that involve a total of 1515 accreted satellites (with 100–150 satellites contributing to each halo) from Bullock & Johnston (2005). Each dark matter host of the 11 halo realizations has a total mass of M_virial(z = 0) = 1.4 × 10$^{12}\;{{M}_{\odot }}$ generated by merger trees using a statistical Monte Carlo method with an extended Press–Schecter formalism (Somerville & Kolatt 1999; Lacey & Cole 1993; Bullock & Johnston 2005, and references therein). Differences in the AHP between each halo are entirely based on the randomness in the merger trees.

CARDs for these 11 merger histories were generated from a semianalytic chemical enrichment code (Robertson et al. 2005) that was applied separately to each infalling satellite and combined with the simulations by Font et al. (2006). Since the enrichment code was implemented for each satellite generated, we can look at individual satellites to assess their relative contribution to their host halo's CARDs. Also, since the aim of this study is to determine the amount of information we can retrieve via chemical abundance observations, we abstain from utilizing any of the satellites' spatial information in our analysis. The main factors contributing to the SFH in the satellites are (1) the epoch of reionization, z_re, (2) the fraction of gas remaining/accreted in the satellite halo after reionization (set mainly by the satellite's virial mass at its time of accretion), (3) the global star formation rate (SFR), and (4) the termination of star formation at the time of accretion (Bullock & Johnston 2005).

Here, one should take note that the assumption of quenched star formation upon satellite accretion into the host galaxy's halo is tentatively supported by observations in Geha et al. (2012) and Grcevich & Putman (2009) (strengthened by Grcevich & Putman 2010). Together, the observations suggest that the star formation of dwarf satellite galaxies is only quenched by close interaction with the host galaxy's halo.⁵

The chemical enrichment of the individual dwarf satellites is affected by these four factors, which are utilized in the simulations to determine the amount of gas available to produce stars and the duration of star formation, which, in turn, determines the chemical evolution of each satellite as prescribed in Robertson et al. (2005). The prescription includes α- and Fe-element enrichment from Type II and Type Ia supernovae (SNe) and stellar wind outflows of metals. The chemical evolution model was tuned with an SN feedback treatment to agree with the local dwarf galaxy stellar mass–metallicity relation (Robertson et al. 2005; see Section 2.3 for further discussion). The α-element patterns in dwarfs versus the smooth halo are consistent with the CARDs of dwarfs found in the compilation of data in Figure 12 of Geisler et al. (2007) (see Figure 1)—an agreement that further bolsters our approach in this investigation (Font et al. 2006).

The function f(x_d) represents the density distribution produced by n random "observations" in chemical abundance space x_d of "stars" (star particles; see Section 2.1 for explanation). Sample distributions for each halo are constructed by randomly drawing "stars" from the halo field.⁶ To mimic observational errors during mock observations, we add a random number drawn from a Gaussian with a dispersion of 0.05 dex to both [α/Fe] and [Fe/H] abundance ratios. The choice of the size of these errors is meant to probe the foreseeable potential of this technique by employing the best-possible conditions for analysis. Evaluation of this technique with ideal conditions provides us with a baseline for expectations from which analysis of real observations in the future can be assessed. In our study, we select samples of $n\approx {{10}^{3}}$, 10⁴, and 3 × 10⁴ representing current, near-future, and optimistically anticipated sample sizes, respectively (K. Freeman 2010, private communication).

Figure 2 shows a 2D CARD ([α/Fe] versus [Fe/H]) of $\sim 3\times {{10}^{4}}$ star particles representing mock stellar abundance ratio observations from the halo 1 simulation. The color of each particle represents the stellar mass of its parent satellite relative to all other accreted satellites. Black and purple particles are donated from the least massive satellites, while orange and red particles are donated from the most massive satellites. The distribution of particles shown demonstrates the expectation that the most massive satellites should account for the vast majority of stars found in the accreted halo stellar population. In comparing this 2D CARD with the observed CARDs in Figure 1, we see that the distributional spread between observed accreted dwarfs of different masses mirrors the distributional spread (in mass) for the simulated dwarfs.

Zoom In Zoom Out Reset image size

The black dashed lines that overlay the colored particle distribution of Figure 2 represent chemical evolution tracks (from the simulations) of typical dwarf masses accreted over the lifetime of the halo. The length of these tracks is primarily affected by the satellite's accretion time. The more time a satellite has to produce stars, the longer its galactic chemical evolution can continue to advance to higher metallicities, and vice versa. The curvature of these tracks is primarily determined by the satellite's mass. The more mass a satellite has to produce stars, the higher its SFR, which means that chemical enrichment by core-collapse SNe is greater. This enhanced early enrichment from core-collapse SNe leads to higher galactic metallicities before the typical 1 Gyr onset (delay) in Type Ia SNe begins (ends), leading to the establishment of a so-called [α/Fe]-knee via significant contributions to Fe abundances. The incorporation of these various tracks into our dwarf model templates is discussed in the next section.

To see whether we can recover the AHP of our simulated halos from our mock observations, we need to generate templates that represent typical accretion events of given satellite stellar mass and age. The most "naive" approach to creating our templates would be to evenly divide the possible range in time t_acc (0–13 Gyr) and mass (stellar) M_sat (${{10}^{0-9}}\;{{M}_{\odot }}$). This division would form N_r mass-binned templates (rows) by N_c time-binned templates (columns) with some "empty" templates (N_empty) where the total number of templates equals ${{N}_{{\rm temps}}}={{N}_{r}}\times {{N}_{c}}-{{N}_{{\rm empty}}}$. However, since decades in galactic (stellar) mass have intuitive implications for galaxy evolution, we restrict our current templates to even divisions in t_acc while we divide M_sat by decades of mass from ${{10}^{5}}$ to ${{10}^{9}}\;{{M}_{\odot }}$ and combine all satellites below ${{10}^{5}}\;{{M}_{\odot }}$ into a fifth mass bin (see Figure 3).

Zoom In Zoom Out Reset image size

After divisions in the t_acc –M_sat plane are selected, all 1515 dwarf satellite models are divided among the bins created by the selected partitions based on each dwarf's individual t_acc and M_sat. During the process, each dwarf's chemical track⁷ (see Section 2.2) is smeared out by a convolution of each star particle with an observational error of ${{\sigma }_{{\rm err}}}=0.05$ dex in both chemical dimensions. To generate the CARDs required for implementation of our recovery algorithm (i.e., the EM algorithm), we separate an average of ∼19,500 star particles per satellite (with errors) into square bins of 0.1 dex that span 3 dex in [Fe/H] (−3 to 0 dex) and 1.7 dex in [α/Fe] (−0.7 to 1). The collection of all binned distributions in our 2D chemical space is normalized to produce an ensemble of probability densities that represent our satellite template set (STS).

Figure 3 shows our 5 × 5 STS as an example of our model template scheme. The full 5 × 5 panel (top right) shows the evenly spaced bins in accretion time versus bins spaced out by decades of accreted satellite stellar mass down to ${{10}^{5}}\;{{M}_{\odot }}$, below which all other satellites are binned together. As stated in Section 2.1, the feedback prescriptions in the chemical evolution models were tuned to reproduce the chemical abundance relationships observed in galactic surveys. First, the mass (luminosity) versus metallicity ([Fe/H]) relationship can be seen by inspecting the trends along any accretion time column. This relationship shows an increase in the distribution peak value of [Fe/H] (and a decrease in the distribution peak value of [α/Fe]) with increasing mass (luminosity) of the galaxy. Second, an age–metallicity relationship can be seen by inspecting the trends along any accreted satellite mass row (i.e., when holding the mass range constant). This relationship shows a decrease in the distribution peak value of [Fe/H] (and an increase in the distribution peak value of [α/Fe]) with an increase in the accretion time epoch (i.e., which dictates the available time for star formation) of the galaxy. However, it should be noted that this age–metallicity relationship is not strictly expected to hold for any given set of dwarf galaxies as other processes are as likely to quench star formation in dwarfs before accretion takes place.

Projections of the 5 × 5 STS, in accreted satellite mass and accretion time, are shown in the top left and bottom right corners of Figure 3, respectively. A comparison of both projections reveals smaller differences in CARDs between adjacent bins of accretion time than in adjacent bins of accreted satellite mass. The similarities between dwarf models in the 5 × 1 STS projection suggest that the EM algorithm will perform better when utilizing the 1 × 5 STS projections of accreted satellite mass for estimates (see Sections 3.2, 3.3, and 4.2 for further discussion). Finally, a 1 × 1 STS projection displaying the probability density function of our master template (i.e., containing the CARDs of all 1515 simulated dwarfs) is shown in the bottom left corner of the figure.

In Section 3 we use the two 1D projections discussed here to form a basis of analysis for the EM algorithm's performance and our ability to recover the AHP of halos in one dimension of mass or time.

The composition of our halos can be best described as a finite mixture of discrete accreted objects that exhibit varying characteristics in a shared CARD space (x = [α/Fe], y = [Fe/H]). Since we can construct models for these accreted objects, we can create a mixture model

$$\begin{eqnarray}&&f({{x}_{i}},{{y}_{i}})=\mathop{\sum }\limits_{j=1}^{m}\;{{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}})\end{eqnarray} \tag{ 2 }$$

where the relations

$$\begin{eqnarray*}&&\mathop{\sum }\limits_{j=1}^{m}{{A}_{j}}=1\ \ \ \ {\rm for}\ {{A}_{j}}\geqslant 0,\ \ j=\{1,...,m\}\end{eqnarray*}$$

confine the relative contribution of model satellites ${\boldsymbol{A}}$. Given n observations of $\{{{x}_{i}},{{y}_{i}}\}$, we can construct a log-likelihood function as follows:

$$\begin{eqnarray}\begin{array}{rcl} L({\boldsymbol{A}} ) & = & \mathop{\prod }\limits_{i=1}^{n}f({{x}_{i}},{{y}_{i}})={{\mathop{\prod }\limits_{i=1}^{n}}^{{}}}\left\{ \mathop{\sum }\limits_{j=1}^{m}{{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}}) \right\} \\ {\rm log} \;L({\boldsymbol{A}} ) & = & \mathop{\sum }\limits_{i=1}^{n}{\rm log} \left( \mathop{\sum }\limits_{j=1}^{m}{{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}}) \right). \\ \end{array}\end{eqnarray} \tag{ 3 }$$

Maximizing ${\rm log} \ L({\boldsymbol{A}} )$ will yield the maximum likelihood estimate ${{{\boldsymbol{A}} }_{{\rm MLE}}}$ for ${{{\boldsymbol{A}} }_{{\rm EM}}}$—our best EM estimate for the true A_j values ${{{\boldsymbol{A}} }_{T}}$. This task, which can be computationally arduous, can be made tractable by adding a latent indicator, z, to each observed data point (x, y), to represent the model template of origin. By designating data set $\{{{x}_{i}},{{y}_{i}},{{z}_{i}}\}_{i=1}^{n}$ as our complete data, we can then define a complete data likelihood as

$$\begin{eqnarray}&&L({\boldsymbol{A}} )={{\mathop{\prod }\limits_{i=1}^{n}}^{{}}}\mathop{\prod }\limits_{j=1}^{m}{{\left\{ {{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}}) \right\}}^{{{z}_{ij}}}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\ell ({\boldsymbol{A}} )=\mathop{\sum }\limits_{i=1}^{n}\mathop{\sum }\limits_{j=1}^{m}\;{{z}_{ij}}\;{\rm log} \{{{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}})\}\end{eqnarray} \tag{ 5 }$$

where z_ij equals the hard expectation that (${{x}_{i}},{{y}_{i}}$) comes from the jth satellite template and $\ell ({\boldsymbol{A}} )$ is the complete data log-likelihood.

As stated above, the log-likelihood derived above can be used to obtain ${{{\boldsymbol{A}} }_{{\rm EM}}}$ via the EM algorithm. Starting from an initial set of guesses, ${{{\boldsymbol{A}} }^{(0)}}$, the algorithm iteratively steps through guesses (which are informed by the former set) until the value of the log-likelihood $\ell ({\boldsymbol{A}} )$, conditioned on the data (and within some tolerance), is maximized. More specifically, the maximizing value of the tth iteration, ${{{\boldsymbol{A}} }^{(t)}}$, is then used as the starting value for the next run, and it continues until the likelihood changes by less than 10⁻³ over 25 iterations.⁸ Details of the implementation of this technique are shown in the Appendix.

We discuss how we evaluate the success of our estimates in the next section. Results from the EM estimates are discussed from Section 3 onward.

In order to evaluate the relative success among our calculated AHPs across all halos and the success of the technique across various STSs, we compare the EM estimates, ${{{\boldsymbol{A}} }_{{\rm EM}}}$, with the known true values, ${{{\boldsymbol{A}} }_{T}}$. Using these values, we can calculate the "factor-of-error" (FoE) ratio for each template EM estimate. The FoE value is defined as the maximum between ${{A}_{{\rm EM},j}}/{{A}_{T,j}}$ and ${{A}_{T,j}}/{{A}_{{\rm EM},j}}$.⁹

One way to evaluate the fidelity of our results is to determine an average FoE ratio ($\langle {\rm FoE}\rangle$) from all FoE measured (i.e., from a given STS and halo). This $\langle {\rm FoE}\rangle$ is an average of all FoE_j, weighted by w_j, and given as

$$\begin{eqnarray}&&\langle {\rm FoE}\rangle =\mathop{\sum }\limits_{j=1}^{m}{{w}_{j}}\cdot {\rm Fo}{{{\rm E}}_{j}}\end{eqnarray} \tag{ 6 }$$

where w_j represents a choice of weights for the relative importance of each template estimate and m is the number of templates used. The lowest $\langle {\rm FoE}\rangle$ value indicates the best results balanced by w_j in STS templates for each halo examined. For our primary analysis we take a mean of FoE values (${{w}_{j}}={{m}^{-1}}$), while other weights are examined in Section 5. The method of evaluation is applied to results in Sections 3–4.1.2.

As discussed in Section 2.3, we can construct a true AHP from our model stellar halos to determine how accurately we can estimate them using the EM algorithm discussed in Section 2.4. Here, we examine the accuracy of our 1 × 5 STS estimates, which are a 1D set of five mass bins (as described in Section 2.3 and shown in the top left corner of Figure 3)—that is to say, we evaluate how well we can recover stellar mass fraction contributions from satellites with no sensitivity to their time of accretion.

Figure 4 presents some characteristic results from our 1 × 5 STS analysis. The top panel indicates that open squares represent the ${{{\boldsymbol{A}} }_{{\rm EM}}}$ values estimated by applying our EM analysis to observed abundances from $\sim {{10}^{4}}$ observed stars.¹⁰ Error bars (calculated from the Fisher information matrix) indicate the smallest possible (1σ) error values (see the Appendix for details). The colored circles shown represent the ${{{\boldsymbol{A}} }_{T}}$ (true) values, while the specific colors of each circle categorize the FoE between ${{{\boldsymbol{A}} }_{{\rm EM}}}$ and ${{{\boldsymbol{A}} }_{T}}$ values by the color legend to the right of the plot. Various FoE values spanning less than 1.1 (bluish-green filled circle) to 10 or more (red filled circle) are examined.

Zoom In Zoom Out Reset image size

In the figure, two plots are chosen to display results from two representative halos (labeled by "h" with the designated number for the halo for short). The two halos are the best (h8) and worst (h5) AHP estimates as determined by their average FoE ($\langle {\rm FoE}\rangle$) values.

Looking at our best EM estimates from h8, we see that individual A_EM produce errors that are within a factor of 2.5 or better for all template estimates using $\sim {{10}^{4}}$ observed stars. This is remarkable considering that we are characterizing $\lesssim {{10}^{-2}}$ to 10⁻³ of the total halo luminosity for the lowest-mass bins.

Our worst EM estimates from h5 seem to reinforce the notion that this analysis provides reliable results. In this worst-case scenario, most estimates are within a factor of two, while the worst estimate (given for our most massive satellite template) is within a factor of eight.

The other principle dimension of our analysis is time. Using the same prescribed analysis above, we can examine the success of estimating AHP from a 1D set of five equally spaced time bins (also described in Section 2.3)—that is to say, we evaluate how well we can recover stellar mass fraction contributions from satellites with no sensitivity to their stellar masses. Figure 5 presents some characteristic results from our 5 × 1 STS analysis. In the figure, plots are chosen based on the same criteria used in making Figure 4. The best EM estimates from h6 reveal very different results concerning the reliability of our analysis when compared to the 1D mass-resolved template results. While both the two most recent and two earliest accretion events have FoE values $\leqslant 2.5$, the "medieval" accretion event has an FoE value $\gtrsim 30$. Here, only the least massive accretion event has a poor FoE value.

Zoom In Zoom Out Reset image size

Our worst EM estimates from h7 follow a trend where all but the most massive accretion event (the medieval event in this case) have markedly poor FoE values that range from $\gtrsim 20$ to ≳10³. Here, the best estimate has an FoE $\leqslant 1.5$ (i.e., with 50% of the true value). While the estimates call into question the reliability of using multiple dimensions in t_acc and M_sat, the overall results were already anticipated from the visual inspection of these templates in the bottom right corner of Figure 3. As suggested earlier, it is likely that degeneracies in CARDs within this template set led to the poor ${{{\boldsymbol{A}} }_{{\rm EM}}}$ estimates seen. In particular, the difference between FoE values for the medieval accretion events in h6 and h7 versus the other events comes down to the dominant accretion event templates subsuming those events that are both highly degenerate in CARD space and significantly less massive than the main event(s). As a consequence, it may appear hopeless to try to glean any information about the accretion times from 1D estimates. This may also hold true for estimates in multiple dimensions when accretion time is treated as the dominant dimension of analysis (see Section 4 for further discussion).

Our complete results, summarized by $\langle {\rm FoE}\rangle$, provide us with insights into the overall effectiveness of the analysis for all 11 halos. Figure 6 displays $\langle {\rm FoE}\rangle$ values for the 1 × 5 STS (i.e., 1D mass-resolved; top panel) and the 5 × 1 STS (i.e., 1D time-resolved; bottom panel). In both panels, each plot shows a histogram of $\langle {\rm FoE}\rangle$ values, calculated using the number of observed stars indicated in each plot, and normalized by the number of halos examined. Dotted light-gray lines indicate an FoE = 2, which indicates, by eye, the vast difference in trying to recover AHPs from 1D mass of accreted satellite templates versus 1D time of accretion templates.

Zoom In Zoom Out Reset image size

In our mass-resolved (1 × 5 STS) estimates (top panel), we can examine the overall success of these templates and note the degree of improvement in estimates as a result of using more data points. Looking at the full panel, we can clearly see the gradual, distinct improvement in ${{{\boldsymbol{A}} }_{{\rm EM}}}$ estimates when a larger data set is used. The median $\langle {\rm FoE}\rangle$ (i.e., our accuracy) for each larger set of observed stars is ∼2.55, ∼2.16, and ∼2.06, respectively. However, it is important to note that the modest improvement between the last two data sets possibly indicates that the method is hitting a limit owing to the number of templates versus the number of stars used. In our time-resolved (5 × 1 STS) estimates (bottom panel), we can see that these estimates are far poorer than the estimates for the mass-resolved estimates. In fact, the time-resolved estimates have a median $\langle {\rm FoE}\rangle$ for each larger set of observed stars equal to ∼100, ∼175, and ∼192, respectively.

Even more critical is the fact that these estimates get marginally worse with number of observations used. This suggests that there are degeneracies between templates in the set that cannot be removed with more CARD information in just two chemical abundance ratio dimensions. Conversely, these degeneracies may also suggest the inherent need for mass divisions in the STS to see differences in templates—a possibility that motivates moving our STS to higher dimensions in the ${{t}_{{\rm acc}}}-{{M}_{{\rm sat}}}$ plane. In the next sections, we discuss the impact of expanding our analysis to multiple dimensions in order to achieve better estimates.

The goal of expanding our STS into higher dimensions is twofold. First, we want to directly recover AHPs with high fidelity by dividing our ${{t}_{{\rm acc}}}-{{M}_{{\rm sat}}}$ plane into templates that would reveal interesting information (e.g., about the MW halo's history) when applied to real abundance observations. Second, we want to indirectly recover 1D stellar mass functions (mass-resolved profiles) and time of accretion histories (time-resolved profiles) of our halos by summing "like" estimates in time or mass together (marginalization) in order to generate better accounts in 1D than could be done directly. Our hypothesis is that allowing a finer grid in time will produce templates with less degeneracy and allow a better recovery of the AHP. Of course, this must be balanced by the size of our sample and its ability to constrain the additional parameters (larger ${{{\boldsymbol{A}} }_{{\rm EM}}}$ set) from the increased number of templates.

4.1.1. "Early" versus "Recent" Accretion: 2 × 5 STS Results

To address our goals, we start by generating templates for our 2 × 5 STS, which have two evenly divided time bins for recent (0–6.5 Gyr ago) and early (6.5–13 Gyr ago) epochs. Figure 7 displays a selection of results that reveal the success of EM estimates due to the application of our 2 × 5 STS. In the figure, we can once again examine the best (h11), the median (h2), and the worst (h7) of the halo estimates using these templates. In the first column of Figure 7 we display the values of $\langle {\rm FoE}\rangle$ to indicate the success of estimates using $\sim {{10}^{4}}$ stars, which can be compared with our marginalized results in the rightmost columns. At first glance, we see that all panels indicate, by (mostly green) colors, that most estimates are within an FoE of 2. For the best EM estimates (from h11), it is encouraging that all FoE values are $\leqslant 2$.

Zoom In Zoom Out Reset image size

However, for the worst EM estimates (from h7) we see a marked decrease in the fidelity of a couple of estimates and especially for one at the high-mass end. Here, we see that the ${{A}_{T,j}}$ value for the early accreted ${{10}^{7-8}}\;{{M}_{\odot }}$ template is actually similar to its recently accreted counterpart, whereas the EM estimates are very different. While the early accretion event is estimated to be essentially nonexistent, both the adjacent higher-mass template (early accreted ${{10}^{8-9}}\;{{M}_{\odot }}$ template) and the recent accretion 10^7–8 $\;{{M}_{\odot }}$ template have slightly higher EM estimates than their true values. The $\leqslant 50\%$ difference in FoE values is probably due to both templates subsuming the contributions from the poorly estimated ${{10}^{8-9}}\;{{M}_{\odot }}$ template. Given that this template is high mass and accreted early, this degeneracy is likely due to the fact that the accretion of most massive systems happens early in most of the 11 halos' histories. Since the 1515 satellites used to make the templates are composed of 11 ensembles of accreted dwarf systems that make up the composition of our simulated halos, it is not surprising that coarse divisions in accretion epochs lead to disparities in the fidelity of our estimates across the 6.5 Gyr divide.

On the other hand, as indicated by our best selection, it is reassuring that, given the simplicity of our dwarf models, there is enough information in their CARDs to make templates that differentiate between higher-mass progenitors of the halo at different epochs. This is true despite the fact that the highest-mass dwarf models show the greatest amount of degeneracy among accreted systems throughout all halos' assembly histories. Also, given the strength of current techniques to more accurately identify recent galaxy formation (e.g., color–magnitude diagrams from photometric surveys that lead to estimates for age and SFHs and phase-space diagrams from low-resolution spectroscopic surveys that lead to estimates for accretion histories), it is encouraging that our technique works so well for early accretion epochs and low-luminosity objects.

In the last column of Figure 7, we present a summation of estimates across accretion epochs (shown with $\langle {\rm FoE}\rangle$ values labeled "L") and across binned satellite luminosities (labeled "T") for all epochs. Here, we confirm that a marginalization of estimates across our two epochs yields 1D estimates with greater fidelity than its 2D decomposition for the worst EM estimates, as indicated by the L-labeled $\langle {\rm FoE}\rangle$ values. More importantly, we can compare our best worst values for our h7 estimates (FoE = 2.161) with the respective 1D h8 estimates (FoE = 2.059) in Figure 4. A comparison of these values shows tentative evidence that our hypothesis about gains in STS information is correct—that the 1D marginalizations across epochs from a 2D STS provide on par or better estimates for 1D AHP than does our bona fide 1D STS. We can also compare the set of "T"-labeled best $\langle {\rm FoE}\rangle$ values for our 1D marginalizations across satellite luminosity bins in Figure 7 with the set of values calculated for Figure 5 (FoE = [7.168, 485.6] for our best and worst values, respectively). Here, we find that our estimates for our time of accretion histories improve substantially overall and dramatically when comparing our best and worst AHP estimates. The next two sections address whether these improvements are ubiquitous as we increase the resolution of our STS in the accretion time dimension.

4.1.2. "Medieval" Accretion: 3 × 5 STS Results

In order to further test our ability to estimate AHPs, we seek to increase our accretion time resolution (by adding an intermediate "medieval" accretion epoch), with the hopes that greater information from an expanded STS will lead to better AHP estimates.

Figure 8 shows our best and worst 3 × 5 STS results. The $\langle {\rm FoE}\rangle$ values between the best and worst EM estimates show a substantial decrease in quality. It is immediately apparent (from color) that individual estimations fared significantly worse than they did in the 2 × 5 STS selections of Figure 7. Also, by inspection, the medieval epoch yields the worst estimates overall. Similar to Figure 7, early epoch estimates of Figure 8 are the most accurate. The overall decrease in performance from our 2 × 5 to 3 × 5 STS is likely due to the degeneracy in CARD space between some adjacent templates in the 3 × 5 STS (e.g., see Figure 3 for illustration of this effect) and across accretion time for the higher-luminosity templates. For example, if we look across the recent and medieval epochs for our worst EM estimate selection, we can see that there are degeneracies in the estimates for the highest stellar mass bins (${{10}^{8-9}}\;{{M}_{\odot }}$). These degeneracies are due to the increasing similarities between chemical model tracks of more massive (and luminous) dwarf satellite models. Such degeneracies can lead to the satisfaction of estimates across all epochs by one individual template (e.g., h7 from Figure 5), by distributing the luminosity fraction among co-degenerate templates (e.g., h7 from Figure 7), or by swapping estimates across adjacent epochs (e.g., h10 from Figure 8). However, it appears that a clear separation in accretion epochs for the same stellar mass bins possibly reduces degeneracies between them (as seen for the best (h3) estimates).

Zoom In Zoom Out Reset image size

If we look at the final column for our 1D marginalizations from the 2D 3 × 5 STS, we once again see improvements in $\langle {\rm FoE}\rangle$ values in comparison to Figures 4 and 5 (e.g., look at "L" and "T" values for all selections vs. uniformly weighed values in Figure 12 in Section 5). While improvements were anticipated, it is still surprising, given the relative lack of success for individual 3 × 5 STS templates, that marginalization of the worst 3 × 5 STS leads to 1D estimates that offer an improvement over the 2 × 5 STS marginalized 1D estimates. In this case, some inaccuracies due to degeneracies across epochs are mitigated by summation over accretion epochs. Consequently, improvements to our marginalized mass-resolved 1D estimates arise from an increase in the STS epoch resolution. Presumably, the better estimates would originate directly from improved individual epoch estimates. However, poor individual estimates due to degeneracies within the same stellar mass bins refute this idea. Indeed, it is more likely that improvements to our epoch resolution led to better estimates indirectly, not by decreasing the degeneracies between adjacent epochs, but rather decreasing degeneracies between adjacent stellar mass bins. While the effects described above are certainly taking place, it is still unclear from Figures 4, 5, 7, and 8 whether these improvements remain across all 11 halos. In the next section we examine the $\langle {\rm FoE}\rangle$ values as ensembles across the 11 halos to determine the overall success of recovering AHPs given our STS.

In this section we compare results from all our simulated halos and the templates we constructed. Using FoE values (see Section 2.5), we can determine a cumulative distribution function (CDF) of FoE values with respect to ${{{\boldsymbol{A}} }_{T}}$ for each STS used. The CDF values described above (which we call ${{A}_{T,{\rm sum}}}$) indicate the fraction of the total stellar halo mass we can identify within a given FoE value.

First, we construct ${{A}_{T,{\rm sum}}}$ values in Figure 9 for 6 of our 10 STSs. Each plot frames the recovery of AHPs in terms of the level of accuracy (i.e., FoE) at which we can characterize a certain portion (${{A}_{T,{\rm sum}}}$) of the total luminous stellar content of the halos examined. Once again, differences in the fidelity of our estimates between 5 × 1 and 1 × 5 STSs are clearly shown with a median ${{A}_{T,{\rm sum}}}$ (fraction recovered) with an FoE ≲ 2 being $\simeq 73\%$ and 95%–99%, respectively. Characterizing the success of the method overall, we find that the median ${{A}_{T,{\rm sum}}}$(with FoE ≲ 2) across most STSs is $\simeq 70\%$ or better. It is evident from the STS shown in Figure 9 that EM estimates fair poorly when applied to certain halo realizations. We discuss possible causes for the often poorer estimates of a few halos in Section 5.

Zoom In Zoom Out Reset image size

Figure 10 displays another way we can summarize our results with the utilization of ${{A}_{T,{\rm sum}}}$ and FoE. In the three panels, box-and-whisker plots illustrate the median and shape of the distribution of ${{A}_{T,{\rm sum}}}$ values calculated for estimates with FoE ≲ 2 among all 11 halos.¹¹ The top panel displays similar information to the results shown in Figure 9. The middle and bottom panels show both genuine and marginalized estimates for the 1 × 5 STS accreted mass functions and the 5 × 1 STS accretion time histories, respectively.

Zoom In Zoom Out Reset image size

In the top panel, ${{A}_{T,{\rm sum}}}$(FoE ≲ 2) is plotted, as a color box, for all STSs examined. Here, as in Figure 6, the color refers to the respective number of observations used (as indicated in the plot legend). In the plot, we see that our best median ${{A}_{T,{\rm sum}}}$ values are given by the 1 × 5 and 2 × 5 STSs, while the worst values are given by 5 × 1 and 7 × 5 STSs. The average among the best and worst ${{A}_{T,{\rm sum}}}$ values across all STSs examined and for an increasing number of stellar observations is ∼0.96–0.98 and ∼0.29–0.41, respectively. The average median ${{A}_{T,{\rm sum}}}$ values across all STSs examined and for an increasing number of stellar observations are 0.742, 0.783, and 0.785. This means that, on average, our FoE are ≲2 for at least ∼75% of the total halo stellar mass (i.e., ${{A}_{T,{\rm max} }}=0.75$) observed.

Marginalized values, which are defined in Section 4.1.1, are useful for evaluating any gains that may potentially arise owing to better time (or mass) resolution. More precisely, any information about templates that is lost or gained should generally result in a corresponding rise or drop in $\langle {\rm FoE}\rangle$ and thus appear as an increase in ${{A}_{T,{\rm sum}}}$(FoE ≲ 2). As a reference, a gray bar is placed in each panel to indicate a region where the ${{A}_{T,{\rm sum}}}$(FoE ≲ 2) values range from 70% to 100% (from bottom to top).

The middle panel shows our mass-resolved marginalized values (summed over accretion time bins) for eight of the nine STSs (with 5 × 1 omitted because its value is not applicable in this context). The plot shows an across-the-board increase in ${{A}_{T,{\rm sum}}}$(FoE ≲ 2) values (i.e., a general drop in all STS $\langle {\rm FoE}\rangle$ values) measured for a recovery of the total stellar mass function. The improvement in FoE values despite the tendency for various individual FoE STS values to increase with an increase in the number of templates used, indicates that significant gains were made by using a larger template set for the specific purpose of generating more accurate estimates of a halo's total stellar mass function (via marginalization).

The bottom panel shows our time-resolved marginalized values (summed over mass bins) for eight of the nine STSs (with 1 × 5 also omitted because its value is not applicable in this context). In this case, the plot shows a descending trend in ${{A}_{T,{\rm sum}}}$(FoE ≲ 2) values with larger STS (i.e., a generally ascending rise in $\langle {\rm FoE}\rangle$ values with increasing STS size) measured for a recovery of the total accretion time history. Despite the decrease in ${{A}_{T,{\rm sum}}}$(FoE ≲ 2) values, these values remain relatively good (above 70% for ${{A}_{T,{\rm sum}}}$ values above the bottom 50% margin) up to our 6 × 5 STS. Indeed, all time-resolved marginalized values show a significant improvement in accretion time histories over the history given by the 5 × 1 STS. Overall, the results show that we could expect to recover accretion time histories using the EM algorithm given that we use reasonable templates.

Results shown in Figures 9 and 10 prove that even with the simplest template divisions, we could, with the appropriate data set, recover the accretion history of the MW halo. To that point, we find that these STS EM estimates can recover the total contributions from accreted systems (templates) of similar mass (i.e., halo luminosity function) to within a factor of 1.02 (≤2% of the true value) for most of the 11 halos. Separately, the EM algorithm can determine the mass fractions within accretion times to within a factor of $\gtrsim 4$ for at least 90% of the halo's total stellar mass. Both results present encouraging prospects for recovering the accretion history of the MW halo from current and near-future data collections.

We can test the statistical robustness of the EM algorithm's application to our simulated halos by performing a likelihood ratio test on the results of our analysis. By determining the true (${{{\boldsymbol{A}} }_{T}}$) and respective ${{{\boldsymbol{A}} }_{{\rm EM}}}$ likelihood values from each application of STS to our halos via the EM algorithm, we can calculate a ${{\chi }^{2}}$-statistic defined by the following equation:

$$\begin{eqnarray}&&{{\chi }^{2}}=-2\ {\rm ln} \left( \frac{{{\lambda }_{T}}}{{{\lambda }_{{\rm EM}}}} \right)\end{eqnarray} \tag{ 7 }$$

where ${{\lambda }_{T}}$ and ${{\lambda }_{{\rm EM}}}$ are the likelihoods for ${{{\boldsymbol{A}} }_{T}}$ and ${{{\boldsymbol{A}} }_{{\rm EM}}}$ values, respectively. One can then reject the assumption that the true AHP templates are well approximated by the STS used if the ${{\chi }^{2}}$-value from Equation (7) is larger than the ${{\chi }^{2}}$-percentile values given k degrees of freedom (k = m_EM - m_T)¹² and a confidence level denoted by α. Figure 11 shows the maximum α-value one can assume for a ${{\chi }^{2}}$-distribution before you have to reject the assumption that suitable AHP templates are chosen. For example, an $\alpha =0.05$ corresponds to a confidence that 95% of all samples taken of a given size are well characterized by the STS in use. Here, we find that out of all sample sizes and STSs used, halos 5, 9, and 10 are by far the worst-characterized halos by our STS divisions. For most STSs used, these halos are ill-matched to the generic STS created in our division scheme and therefore challenge the robustness of this method. Such challenges need to be address before this method can be utilized to model the AHP of the MW halo. The solution resides in the development and incorporation of sufficiently realistic models of dwarf CARDs into this method—a goal that will be addressed in future work.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another consideration in assessing the reliability of our method is to determine how well it uncovers AHPs based on the satellite mass regime we are interested in. Taking Equation (6) from Section 2.5, we can calculate $\langle {\rm FoE}\rangle$ values with different weights—i.e., uniform (mean), low-mass preferred, or high-mass preferred—based on what satellite population(s) one prefers to recover. Figure 12 shows the median $\langle {\rm FoE}\rangle$ among all halos for each STS used. The same colors from Figure 10 are used to indicate the number of stars used for the analysis, and symbols and corresponding lines refer to the type of weighting used (see figure legend). Uniformly weighted $\langle {\rm FoE}\rangle$ values are weighted by ${{{\rm m}}^{-1}}$ (i.e., by the number of templates used) and identical to the weighting used for the main results of this paper. Weights that emphasize more accuracy in low- or high-mass satellite AHPs are weighted by the corresponding upper bin mass limits and their reciprocals, respectively.

In the figure, we can see that $\langle {\rm FoE}\rangle$ values for low-mass satellite recovery fair the best, whereas uniform and high-mass satellite recovery-emphasized weights are a factor of $\gtrsim 10$ in all but the three smallest template sets. In other words, when one emphasizes the accurate recovery of low-mass satellites, the weighting favors templates with lower FoE values, which yields lower overall $\langle {\rm FoE}\rangle$ values. This result further clarifies the immediate strengths of the method: it is adept at differentiating between accreted dwarfs of low mass in CARD space owing to the lack of degeneracies in their occupied region of space. Meanwhile, it is clear that while degeneracies exist in the CARD space occupied by high-mass satellites and larger STSs, we are encouraged by the fact that the introduction of more templates can significantly decrease degeneracies in only two dimensions of CARD space.

It is clear from both our results and our reliability tests that the current method fails often for 3 of the 11 halo simulations. From our examination of these three problematic halos we find that all of them show predominately early accretion of massive dwarf galaxies with integrated CARDs that appear to be highly degenerate when compared with the other eight halos' AHP CARDs examined.

To address the degeneracies that exist (particularly among high-mass systems), we posit that differences between mass-dependent (nucleosynthetic) yields for different nucleosynthetic sites and element groups (e.g., see Lee et al. 2013) can be exploited to greatly reduce or remove such degeneracies by expanding the CARD-space basis set.

For example, we only looked at two dimensions in CARD space, whereas more recent work on "chemical tagging" expands the number of dimensions available by establishing the best chemical abundance signatures to pursue in chemical abundance space in order to optimize survey efforts (e.g., the GALAH survey). One way to optimize our surveys for searches in chemical abundance space is to prioritize spectroscopic observations for elements that confer the greatest amount of distinction between systems with different origins. To this end, principle component analysis was used by Ting et al. (2012) to identity and rank the six to nine most distinguishing elements in chemical abundance space. In their work, the chemical abundance space of various parts of both the galaxy and the galactic neighborhood was examined to determine the best elements to observe in order to decipher their galactic chemical evolution. A CARD-space basis set derived from various combinations of these elements is likely to offer the breaks in degeneracies that we require. Lastly, it should be noted that the current and upcoming surveys that are best poised to provide the data required for our approach are the Subaru prime focus spectrograph (PFS; Takada et al. 2014) and the Gaia-ESO (Gilmore et al. 2012) surveys.

To implement the algorithm, we first need to derive the expression for the complete data log-likelihood, given by Equation (5), which is conditioned on the data. To do this, it is necessary to decide on a mode of usage for z_ij. The use of z casts the EM algorithm as either hard when its value discretely indicates the ${{f}_{j}}({{x}_{i}},{{y}_{i}})$ of origin or soft when its value probabilistically indicates the origin of point (${{x}_{i}},{{y}_{i}}$) across all f_j. For this application, we chose to implement a hard EM algorithm for estimation of A_MLE in which z_ij has a true value equal to 1 if the data point (${{x}_{i}},{{y}_{i}}$) comes from model f_j and 0 otherwise. Thus, our overall expectation is

$$\begin{eqnarray}\begin{array}{rcl} {{E}_{{\boldsymbol{A}} }}\left[ \ell ({\boldsymbol{A}} )|{\boldsymbol{x}} ,{\boldsymbol{y}} \right] & = & \mathop{\sum }\limits_{i=1}^{n}\mathop{\sum }\limits_{j=1}^{m}{{E}_{{\boldsymbol{A}} }}[{{z}_{ij}}|{{x}_{i}},{{y}_{i}}] \\ {} & {} & \times \{{\rm log} {{A}_{j}}+{\rm log} {{f}_{j}}({{x}_{i}},{{y}_{i}})\} \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

where

$$\begin{eqnarray}&&{{E}_{{\boldsymbol{A}} }}[{{z}_{ij}}|{{x}_{i}},{{y}_{i}}]=\frac{{{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}})}{\mathop{\sum }\limits_{k=1}^{m}{{A}_{k}}{{f}_{k}}({{x}_{i}},{{y}_{i}})}\end{eqnarray} \tag{ A.2 }$$

as defined by Equation (2). Since we are ultimately maximizing Equation (A.1), the nonconstant term, Equation (A.2), becomes the component of interest. To iteratively evaluate this expectation, we let $w_{ij}^{(t)}$ be Equation (A.2) at the tth step:

$$\begin{eqnarray*}w_{ij}^{(t+1)}=\left\{ \begin{array}{ccccccccccccccc} \frac{{{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}})}{\sum _{k=1}^{m}{{A}_{k}}{{f}_{k}}({{x}_{i}},{{y}_{i}})} & j=1,\cdots ,m \\ 1-{{w}_{i1}}-...-{{w}_{i,m-1}} & j=m. \\ \end{array} \right.\end{eqnarray*}$$

Since ${\boldsymbol{A}}$ is not defined for the first evaluation, we use a random initialization to generate ${\boldsymbol{w}} _{j}^{(0)}$. Here, it should be noted that convergence is not sensitive to the choice of values in our case, though it can be in cases where the likelihood is riddled with local maxima. If we examine the expression above, we can conceptually define the mechanism for maximization as a "ratcheting up" of ${{E}_{A}}[{{z}_{ij}}|{{x}_{i}},{{y}_{i}}]$ values by maximizing ${{A}_{j}}{{f}_{j}}({{x}_{i}},{{y}_{i}})$ with respect to $\sum _{k=1}^{m}{{A}_{k}}{{f}_{k}}({{x}_{i}},{{y}_{i}})$. Derivation of the maximization expression is discussed below.

Above we defined an explicit formulation for the expected log-likelihood (Equation (A.2)) given a single parameter ${\boldsymbol{A}}$ and the data (${\boldsymbol{x}} ,\;{\boldsymbol{y}}$). The argument of the maximum of Equation (A.2) at each iteration t provides an estimate that approaches the MLE of ${\boldsymbol{A}}$ and is given by

$$\begin{eqnarray}&&{{{\boldsymbol{A}} }^{(t)}}=\mathop{{\rm argmax}}\limits_{{\boldsymbol{A}} }\left[ \ell ({\boldsymbol{A}} )|{\boldsymbol{x}} ,{\boldsymbol{y}} ,{{{\boldsymbol{A}} }^{(t-1)}} \right].\end{eqnarray} \tag{ A.3 }$$

Accounting for the m – 1 free parameters of ${\boldsymbol{A}}$, differentiation of Equation (A.1) with Equation (A.2) proceeds, for $k=1,\cdots ,m-1$, as

$$\begin{eqnarray*}\begin{array}{rcl} \frac{\partial }{\partial {{A}_{k}}}{{E}_{A}}\left[ \ell ({\boldsymbol{A}} )|{\boldsymbol{x}} ,{\boldsymbol{y}} \right] & = & \mathop{\sum }\limits_{i=1}^{n}\left\{ w_{ik}^{(t-1)}\frac{1}{{{A}_{k}}} \right. \\ {} & {} & -\left. w_{im}^{(t-1)}\frac{1}{1-{{A}_{1}}-...-{{A}_{m-1}}} \right\} \\ \end{array}\end{eqnarray*}$$

where the first term in the summation accounts for all values of $k\leqslant m$ and the second term eliminates overcounting of the first term at k = m. The derivative of an argmax is always equal to zero since we are taking a derivative at the maximum point of the function in question (in our case the expectation of the log-likelihood). Thus, we can expand the summation of data points and equate the terms described above to one another:

$$\begin{eqnarray*}&&\frac{1}{{{A}_{k}}}\mathop{\sum }\limits_{i=1}^{n}w_{ik}^{(t-1)}=\frac{1}{1-{{A}_{1}}-...-{{A}_{m-1}}}\mathop{\sum }\limits_{i=1}^{n}w_{im}^{(t-1)}.\end{eqnarray*}$$

Consequently, these terms being equal means that every $k\leqslant m$ term is equal to each as shown below,

$$\begin{eqnarray*}&&\frac{1}{{{A}_{k}}}\mathop{\sum }\limits_{i=1}^{n}w_{ik}^{(t-1)}=...=\frac{1}{{{A}_{m-1}}}\mathop{\sum }\limits_{i=1}^{n}w_{i,m-1}^{(t-1)}=c\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\boldsymbol{A}} _{k}^{(t)}=\frac{\mathop{\sum }\limits_{i=1}^{n}w_{ik}^{(t-1)}}{c}\end{eqnarray*}$$

where c is some constant.

The unknown constant c appears problematic, but because $\sum _{j=1}^{m}{{A}_{j}}=1$, algebraic manipulation reveals that c = n, yielding a final solution that can be numerically evaluated:

$$\begin{eqnarray}&&{\boldsymbol{A}} _{k}^{(t)}=\frac{\mathop{\sum }\limits_{i=1}^{n}w_{ik}^{(t-1)}}{n}\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{\boldsymbol{A}} _{m}^{(t)}=1-{{A}_{1}}-...-{{A}_{m-1}}.\end{eqnarray} \tag{ A.5 }$$

Finally, to implement this algorithm, we simply compute an initial value for ${\boldsymbol{A}}$, inserting each component, ${{A}_{j}}$, into a w_ik^t equal to Equation (A.2) (i.e., with k initially identical to j) and then compute that expression with Equation (A.4) to calculate each new corresponding ${{A}_{k}}$. This process is repeated until our iteration criterion is met.

In our case, computation of ${\boldsymbol{A}} \to {{{\boldsymbol{A}} }_{{\rm EM}}}$ converges relatively quickly for all starting values: on the order of 600 iterations, or half a minute, for n = 1000 (given our stopping criteria). Large ${{A}_{{\rm EM},k}}$ values typically emerge after two or three iterations, and most change, absolutely speaking, occurs in the first 50 to 100 iterations. For error estimation, we can provide values for the minimum error possible through an inversion of the Fisher information matrix (see Section A.3 for brief derivation). Although we have an idea of what the best possible errors are, such values exclude the use of more standard approaches to assessments of parameter estimation, like the reduced ${{\chi }^{2}}$ statistic.

The asymptotic covariance matrix of ${{{\boldsymbol{\hat{A}}} }_{{\rm EM}}}$ can be approximated by the inverse of the observed Fisher information matrix, I.

As ${{A}_{{\rm EM},m}}=1-\sum _{j=1}^{(m-1)}{{A}_{{\rm EM},j}}$, there are only $m-1$ free parameters. Thus, let ${\boldsymbol{A}} _{{\rm EM}}^{\prime }=({{A}_{{\rm EM},1}},\ldots ,{{A}_{{\rm EM},(m-1)}})$. Using ${{f}_{ij}}={{f}_{j}}({{x}_{i}},{{y}_{i}})$ for brevity, the likelihood can then be expressed as

$$\begin{eqnarray}\begin{array}{rcl} \ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime }) & = & \mathop{\sum }\limits_{i=1}^{n}{\rm log} \left\{ \left( \mathop{\sum }\limits_{j=1}^{m-1}{{A}_{{\rm EM},j}}{{f}_{ij}} \right) \right. \\ {} & {} & +\left. (1-{{A}_{{\rm EM},1}},\ldots ,{{A}_{{\rm EM},(m-1)}}){{f}_{im}} \right\} \\ \end{array}\end{eqnarray} \tag{ A.6 }$$

The observed information matrix, I, is the $(m-1)\times (m-1)$ negative Hessian of Equation (A.6), evaluated at the observed data points:

$$\begin{eqnarray*}\begin{array}{rcl} I({\boldsymbol{A}} _{{\rm EM}}^{\prime }|{\boldsymbol{x}} ,{\boldsymbol{y}} ) & = & -\frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {\boldsymbol{A}} _{{\rm EM}}^{\prime }\partial {\boldsymbol{A}} _{{\rm EM}}^{\prime T}} \\ {} & = & -\left[ \begin{array}{ccccccccccccccc} \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{{{\partial }^{2}}{{A}_{{\rm EM},1}}} & \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {{A}_{1}}\partial {{A}_{2}}} & \ldots & \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {{A}_{1}}\partial {{A}_{(m-1)}}} \\ \vdots & \vdots & {} & \vdots \\ \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {{A}_{(m-1)}}\partial {{A}_{1}}} & \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {{A}_{(m-1)}}\partial {{A}_{2}}} & \ldots & \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{{{\partial }^{2}}{{A}_{{\rm EM},(m-1)}}} \\ \end{array} \right] \\ \end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl} \frac{\partial \ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {{A}_{{\rm EM},k}}} & = & \mathop{\sum }\limits_{i=1}^{n}\frac{{{f}_{ik}}-{{f}_{im}}}{\mathop{\sum }\limits_{j=1}^{m}{{A}_{{\rm EM},j}}{{f}_{ij}}}\quad {\rm and} \\ \frac{{{\partial }^{2}}\ell ({\boldsymbol{A}} _{{\rm EM}}^{\prime })}{\partial {{A}_{{\rm EM},k}}\partial {{A}_{{\rm EM},r}}} & = & -\mathop{\sum }\limits_{i=1}^{n}\frac{({{f}_{ik}}-{{f}_{im}})({{f}_{ir}}-{{f}_{im}})}{{{(\mathop{\sum }\limits_{j=1}^{m}{{A}_{{\rm EM},j}}{{f}_{ij}})}^{2}}} \\ \end{array}\end{eqnarray*}$$

with $1\leqslant r\leqslant m-1$ such that (k, r) represents the index of the observed information matrix I.

The observed information matrix of ${\boldsymbol{A}} _{{\rm EM}}^{\prime }$ yields the following estimates for covariance and correlation for all m estimated weights in ${{{\boldsymbol{\hat{A}}} }_{{\rm EM}}}$:

$$\begin{eqnarray*}{\rm Cov}({{\hat{A}}_{{\rm EM},p}},{{\hat{A}}_{{\rm EM},q}})=\left\{ \begin{array}{ccccccccccccccc} {{\left[ {{I}^{-1}}({\boldsymbol{\hat{A}}} _{{\rm EM}}^{\prime }) \right]}_{Eq}} & p,q\lt m \\ -\mathop{\sum }\limits_{j=1}^{m-1}{\rm Cov}({{{\hat{A}}}_{{\rm EM},j}},{{{\hat{A}}}_{{\rm EM},q}}) & p=m,q\lt m \\ \mathop{\sum }\limits_{j=1}^{m-1}\mathop{\sum }\limits_{k=1}^{m-1}{\rm Cov}({{{\hat{A}}}_{E,j}},{{{\hat{A}}}_{{\rm EM},q}}) & p,q=m \\ \end{array} \right.\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\rm Var}({{\hat{A}}_{{\rm EM},j}})=\sigma _{j}^{2}={{\left\{ {\rm Cov}({{{{\boldsymbol{\hat{A}}} }}_{{\rm EM}}}) \right\}}_{jj}}\end{eqnarray*}$$