\section{Introduction}
Star clusters are important objects for understanding the formation and evolution of their host galaxies. It is considered that most of star formation is clustered \citep{Lada2003}, therefore knowledge of the physical parameters of a cluster population (e.g., age, mass, chemical composition, and extinction) is essential for constraining the star formation history of the galaxy.

The commonly used method for deriving the physical parameters of unresolved star clusters is based on comparing of their integrated broad-band photometry colors to the colors of simple stellar population (SSP) models; see, e.g., \cite{Anders2004} and \cite{Bridzius2008}. However, this method is strongly biased by the presence of two major problems:
\begin{itemize}
\item degeneracy between various physical parameters of star clusters, see, e.g., \cite{Worthey1994}, \cite{Bridzius2008}. For example, a young cluster possessing high extinction can have colors similar to an older object without extinction -- an age-extinction degeneracy;
\item stochasticity, which is due to the random presence of a few bright stars, which dominate the integrated photometry of unresolved clusters; see, e.g., \cite{Santos1997}, \cite{Deveikis2008}, \cite{MaizApellaniz2009}. Young clusters can have supergiant stars that significantly redden their integrated colors. By using the SSP method, a much older age would be determined for these clusters.
\end{itemize}

Although \cite{Cervino2006} have attempted to describe the problem of stochasticity analytically, the nowadays preferred approach (e.g., \cite{Popescu2009, Popescu2010}; \cite{Fouesneau2010,Fouesneau2012}; \cite{Asad2012}) is to use a Monte-Carlo sampling of stellar IMF to model integrated colors of clusters and build an extensive grid of models, to cover all possible age, mass and extinction ranges. Physical parameters of star clusters are then derived by comparing observations to the grid of models.

Recently \cite{Asad2012} have derived the age and extinction of star clusters in the Large Magellanic Cloud (LMC) using the method of \citet{Popescu2010} that takes stochasticity into account, and compared results to previous derivations based on an isochrone fit to the resolved color-magnitude diagrams (CMDs). They managed to constrain the age of clusters similar to the values found by the isochrone fit only when previously known extinctions of individual clusters were used.

In this paper, a method of deriving physical parameters is developed and applied to star cluster samples (integrated multi-band photometry) in two Local Group galaxies: 1) M31, using a catalog by \cite{Vansevicius2009}, who derived cluster parameters using SSP method, and 2) M33, using objects common to catalogs of \cite{Ma2012} for photometry and \cite{SanRoman2009}, who observed in resolved conditions with the {\it Hubble Space Telescope} ({\it HST}) and estimated age, mass, and extinction based on an isochrone fit to cluster CMDs.

Using this data and artificial cluster simulations, we demonstrate that when the intrinsic range of extinction within the host galaxy is rather narrow, it is not necessary to know the exact value of the extinction for individual clusters to derive their physical parameters reliably.

The paper is organized into the following sections. Section \ref{sec:method} introduces our method of deriving physical parameters of clusters when stochasticity is taken into account, \S\,\ref{sec:grid} describes a grid of simulated cluster models, \S\,\ref{sec:artificial} presents tests of the method on artificial cluster samples, and \S\,\ref{sec:application} applies the method to the M31 and M33 star clusters.

\section{Method of deriving age, mass, and extinction of star cluster}
\label{sec:method}
There are presently two main methods used to derive physical parameters of star clusters based on a 3D grid (age, mass and extinction) of models, which take stochasticity effects into account. The first is a fast $\chi^{2}$ minimization approach used by, e.g., \cite{Popescu2010} and \cite{Beerman2012}, the second, a more accurate but much slower approach, which builds probability maps in the age-mass-extinction space by exploring all the nodes of the grid and selects the most probable solution    \citep[see e.g., ][]{Fouesneau2010}.

The method presented here also explores parameter probability maps; however, by restricting the analysis to the models located in the vicinity of the observed absolute magnitudes (the distance to the object has to be known), we significantly reduce the computation time. The scheme of the method is sketched in Fig.\,\ref{fig1}.

\begin{figure}
\centering
\includegraphics[width=88mm]{fig1.pdf}
\caption{\small Scheme of the method for deriving physical parameters of star cluster. Panel (a) shows the grid of models (only age and mass are displayed; same for extinction). These models are represented in the photometric space (only $M_{U}$ and $M_{B}$ are shown) displayed in panel (b; only 100 models per node without extinction are shown). In panel (c) all the models situated within 3--$\sigma$ (circle) from the observed absolute magnitudes are selected. Panel (d) shows the nodes associated with these models. The 1D age and mass distributions of the selected models are displayed as normalized histograms in panels (e) and (f), used to determine most probable values, indicated by solid vertical lines; confidence intervals are indicated by dashed lines.}
\label{fig1}
\end{figure}

A 3D grid of cluster models is built for every value of the three physical parameters $t,\,m,\,E(B-V)$\,\footnote{We refer to extinction as $E(B-V)$ hereafter.}; for simplicity Fig.\,\ref{fig1}\,(a) shows only a grid for age and mass. For the description of the grid, see \S\,\ref{sec:grid}. Each node of the grid contains 1\,000 models of the same age, mass, and extinction. They populate the photometric parameter space (absolute $UBVRI$ magnitudes). Figure \ref{fig1}\,(b) shows $M_{U}$ and $M_{B}$ with only 100 models per node without extinction, and the used model grid is much more continuous in photometric parameter space.

When the observations are considered in Fig.\,\ref{fig1}\,(c), along with their error bars ($\sigma$; hereafter we use $\sigma = 0.05$ mag for all the passbands of artificial and real cluster samples studied in this paper), which in general can be different for every magnitude, all the models situated within 3--$\sigma$ from the observed magnitudes are selected. Fig.\,\ref{fig1}\,(d) shows the nodes to which the selected models are associated. Other nodes do not play any role in the derivation of parameters, therefore the speed of the algorithm is increased significantly. Finally, the distributions of age and mass, displayed in Figs. \ref{fig1}\,(e, f), are derived from the selected models, as for the extinction, which is not shown in the figure.

For the selected models (Fig.\,\ref{fig1}\,(c) circle) we apply weights as follows: for the models located within 1--$\sigma$ from the observed magnitude a weight of 0.68 is assigned, for the ones between 1--$\sigma$ and 2--$\sigma$ -- 0.28, and for the ones between 2--$\sigma$ and 3--$\sigma$ -- 0.04. The probability density distributions displayed in Figs. \ref{fig1}\,(e, f) are derived by normalizing the total area of each histogram to 1. The solution is taken as the maximum of these 1D distributions. We compute confidence intervals (error bars) by excluding the first and the last 16\% of the area in histograms, following the method of ``central interval'' presented in \S\,2.5.1 of \cite{Andrae2010}.

\section{The age-mass-extinction grid of models}
\label{sec:grid}
To derive physical parameters of star clusters with the method described in \S\,\ref{sec:method}, a large age-mass-extinction grid of models is computed, by applying the algorithm described by \cite{Deveikis2008}. The stellar masses are generated randomly sampling the IMF \citep[corrected for binaries; ][]{Kroupa2001} and their luminosities are derived from stellar isochrone of the selected age and metallicity of the cluster model. The process is continued until the total mass of generated stars reaches the mass of the cluster model. Then, taking the distance to the cluster into account, the magnitudes are computed using the Johnson-Cousins $UBVRI$ photometric system \citep{MaizApellaniz2006}.

For stellar models, we took the PADOVA isochrones\footnote{PADOVA isochrones from ``CMD 2.4'': http://stev.oapd.inaf.it/cmd} from \cite{Marigo2008}, with corrections by \cite{Girardi2010} for the TP-AGB phases. The model grid for a single metallicity $Z=0.008$ contains the following nodes: ages from $\log(t/{\rm yr})$ = 6.6 to 10.1 in steps of 0.05, masses from $\log(m/\mbox{${\rm M}_{\sun}$})$ = 2 to 5 in steps of 0.05. This gives 71 values of age and 61 values of mass, and the grid consists of 1\,000 models per node, i.e. $\sim$ $4\times10^{6}$ stochastic models. To limit the number of models to store in computer's memory, extinction is computed when the observed cluster is compared with the grid of models. It ranges from $E(B-V)$ = 0 to 1 in steps of 0.02, therefore 51 values for the extinction.

To accelerate computation of integrated magnitudes of stochastic star cluster models, we defined a threshold in the isochrone, under which the total luminosity of fainter stars is computed by integration of the stellar luminosity function along the isochrone. Above the threshold, which is defined to select 20\% of the most massive stars, the contribution of the high-mass stars is modeled by the algorithm of \cite{Deveikis2008}. The models built by this improved procedure share the same properties as models built by simulating all stars of the cluster, but require much less computation with a speed gain of a factor $\sim$10.

\section{Test of the method on an artificial cluster sample}
\label{sec:artificial}
\subsection{Degeneracies in an artificial cluster sample}
We simulated artificial star cluster samples with known age, mass, and extinction and used them as input clusters to evaluate the ability of our method to derive physical parameters. The artificial cluster samples consist of 10\,000 clusters with a random age in the $\log(t/{\rm yr})$ range of [6.6, 10.1]. To simulate the mass of input clusters, we used a power-law cluster mass function with index $-2$ in the range of $\log(m/\mbox{${\rm M}_{\sun}$})$ = [2.7, 4.3], so as to have more low-mass clusters in the sample. We have two artificial samples: one without extinction and the other with $E(B-V)$ in the range of [0.0, 1.0] using the Milky Way standard extinction law from \cite{Cardelli1989}.

Figure \ref{fig2} displays the results for artificial cluster sample built without extinction. Panels (a) and (b) show results of the age and mass of the artificial cluster sample without introducing photometric observation errors. Panels (c) and (d) show the case with Gaussian photometric errors of 0.05 mag randomly added to each magnitude of the sample clusters. The introduction of photometric errors results in broadening around the one-to-one line.

In Fig.\,\ref{fig2}\,(b), the asymmetry observed in the mass derivation slightly favors high masses. This is because that in the grid of star cluster models, the nodes of models with higher mass have magnitudes that are less dispersed than the nodes with lower mass models, as a consequence of stochasticity \cite[see e.g., ][]{Deveikis2008}. Thus, when the models of two nodes of different masses are located in the $UBVRI$ ``sphere'' around the observation shown in Fig.\,\ref{fig1}\,(c), the more massive node, with less dispersed models, will dominate. To balance out the effect, a cluster mass function could be used to decrease the importance of the nodes of massive cluster models. Although Fig.\,\ref{fig2}\,(b) shows asymmetry, the means and standard deviations given at $\log(m/\mbox{${\rm M}_{\sun}$})$ = 3.0, 3.5, and 4.0 indicate that this phenomenon is slight.

\begin{figure}
\centering
\includegraphics[width=88mm]{fig2.pdf}
\caption{\small Derived parameters of 10\,000 artificial clusters without extinction. Panels (a) and (b) show age and mass for a sample without photometric errors; (c) and (d) for a sample with Gaussian photometric errors of 0.05 mag. Panel (b) also shows the means and standard deviations (circles with error bars) of the derived mass distributions for clusters with true mass $\log(m/\mbox{${\rm M}_{\sun}$})$ = 3.0, 3.5 and 4.0. The density scale is logarithmic.}
\label{fig2}
\end{figure}

Figure \ref{fig3} gives the results of artificial clusters in the case of extinction. Panels (a), (b), and (c) show the age, mass and extinction derived when there is no photometric errors, and panels (d), (e), and (f) when there are Gaussian photometric errors of 0.05 mag included. In Fig.\,\ref{fig3}\,(d), in the case of photometric errors and extinction, we observe that a broadening around the one-to-one relation increases, especially for $\log(t/{\rm yr}) \gtrsim 8$, which is associated with broadening in the extinction (panel f), which is a sign of the age-extinction degeneracy. It creates two streaks above and below the one-to-one relation in the range of $8 \lesssim \log(t/{\rm yr}) \lesssim 9.5$ that were already perceptible in the case without photometric errors in Fig.\,\ref{fig3}\,(a). However, including of photometric errors does not significantly affect the derivation of mass (panels b and e). We note that a gap in derived ages at $\log(t/{\rm yr}) = 9.15$ is a feature of isochrone due to the increase in the production rate of AGB stars, which was discussed in \cite{Girardi1998}.

\begin{figure*}
\centering
\includegraphics[width=180mm]{fig3.pdf}
\caption{\small Derived parameters of 10\,000 artificial clusters with extinction randomly chosen in the range of $E(B-V)$ = [0, 1]. Panels (a), (b) and (c) show age, mass, and extinction for a sample without photometric errors; (d), (e), and (f) for a sample with Gaussian photometric errors of 0.05 mag.}
\label{fig3}
\end{figure*}

These degeneracy streaks have first been reported by \cite{Fouesneau2010}. The degeneracy streak above the one-to-one line seen in Fig.\,\ref{fig3}\,(d) concerns clusters that are young and that possess intrinsically high extinction, but are derived by the method as older and having lower extinction. Conversely, the streak of clusters below the one-to-one relation involves objects that are derived as older and that have lower extinction than they do in reality. These features are important to keep in mind when deriving the physical parameters of unresolved star clusters.

The treatment of the metallicity effects on the derivation of age, mass, and extinction of star clusters is the subject of a forthcoming paper. To illustrate the effect, we derived the parameters of a cluster sample of $Z = 0.008$ metallicity, successively with a model grid of much lower metallicity, $Z = 0.00013$, and another one with much higher metallicity, $Z = 0.03$. In the first case, the ages are systematically overestimated by $\sim$0.5 dex, the upper streak is more developed, and the lower one disappears. The masses are also overestimated by $\sim$0.5 dex, and there is a preference for overestimated extinction. In the second case, the ages are only slightly ($\sim$0.1 dex) underestimated. The upper streak decreases without vanishing, and the lower one becomes more populated. The masses are underestimated of $\sim$0.1 dex, and the extinction is almost unaffected.

\subsection{Is it possible to reduce the age-extinction degeneracy?}
The upper and lower streaks in Fig.\,\ref{fig3}\,(d) suggest that if a wide extinction range is allowed in a simulated sample, then there are possibilities that a cluster mimics an older one with lower extinction, or inversely a younger one with higher extinction. If the true extinction range of the cluster population is narrow, then we could restrict the search for the extinction within a narrow range in the model grid, resulting in decrease of age-extinction degeneracy.

In Fig.\,\ref{fig4} we show results of the tests for a sample of 10\,000 artificial clusters with true extinction from a range of $E(B-V)$ = [0, 0.5] and with Gaussian photometric errors of 0.05 mag randomly added to each magnitude of the sample clusters. This cluster sample was studied twice, with different allowed extinction ranges in the model grid.

\begin{figure*}
\centering
\includegraphics[width=180mm]{fig4.pdf}
\caption{\small Derived parameters of 10\,000 artificial clusters with true extinction randomly chosen in the range of $E(B-V)$ = [0, 0.5] with Gaussian photometric errors of 0.05 mag. Panels (a), (b), and (c) show age, mass and extinction when a sample is studied using a grid of models with an allowed extinction in the range of [0, 1]; (d), (e), and (f) show results for the same cluster sample studied with a narrower allowed extinction in the range of [0, 0.5]. Panels (c) and (f) are denser than in Fig.\,\ref{fig3}, because the same number of clusters are located on a smaller range of extinction. A square indicated in panels (a) and (d) is used to quantify the degeneracy effect.}
\label{fig4}
\end{figure*}

In the first test, the extinction of the model grid was allowed to vary in a wide range, $E(B-V)$ = [0, 1], shown in Fig.\,\ref{fig4}\,(a, b, c). If a cluster has a true extinction of 0.5, then the maximum underestimation of its extinction can be from 0 to 0.5. But if a cluster has true extinction of 0, the maximum overestimation of its extinction could range from 0 to 1. This explains why the lower streak (i.e. clusters with overestimated extinction) is more extended than the upper one in panel (a).

In the second test, the allowed extinction range of the model grid was reduced to a range of [0, 0.5]; Fig.\,\ref{fig4}\,(d, e, f). The constraint on the extinction range resulted in a reduction of the lower degeneracy streak, seen in Fig.\,\ref{fig4}\,(d), which is less developed than in Fig.\,\ref{fig4}\,(a). From the comparison of Figs.\,\ref{fig4}\,(b) and (e), we see that the mass is less affected by the degeneracy. We note that only the lower streak is modified, since only the higher limit of the allowed extinction range was changed from 1.0 to 0.5. The upper streak is not modified, because the lower limit of the allowed extinction range was not changed.

To quantify the reduction of the lower degeneracy streak due to reduction of the extinction range, all the models situated in the square shown in Figs.\,\ref{fig4}\,(a) and (d) were counted. There are $\sim$480 clusters in that region for the wide extinction range (panel a) and only $\sim$110 for the low extinction range (panel d). The reduction of the extinction range from $E(B-V)$ = [0, 1] to [0, 0.5] thus decreases more than four times the strength of the degeneracy streak.

We conclude that to reduce the degeneracy streaks seen in Fig.\,\ref{fig3}\,(d), we should make a reasonable assumption on the extent of the possible extinction range within the galaxy hosting the studied cluster population.

\section{Application of the method to real star clusters}
\label{sec:application}
\subsection{The M31 galaxy star cluster sample}
A star cluster catalog of 285 objects located in the south west field of the M31 galaxy was compiled by \cite{Narbutis2008} using the deep $BVRI$ and $H_{\alpha}$ photometry from the Subaru telescope, as well as multiband maps based on {\it HST}, {\it Galex}, {\it Spitzer}, and 2MASS imaging. The $UBVRI$ photometry was derived using the Local Group Galaxy Survey data from \cite{Massey2006}. The magnitude limit of the cluster sample was set to the $V = 20.5$ mag.

\cite{Vansevicius2009} selected 238 clusters from that sample, excluding the ones with strong $H_{\alpha}$ emission, and compared their multiband colors to PEGASE \citep{Fioc1997} SSP models to derive their age, mass, metallicity, and extinction. {\it Spitzer} data was used to constrain the maximum extinction for each cluster. They reported $\sim$30 classical globular clusters with low metallicity older than 3 Gyr. The remaining $\sim$210 younger clusters were classified as objects belonging to the disk, with average metallicity $Z=0.008$. Figure \ref{fig5} shows the $U-B$ vs $B-V$ and $U-V$ vs $R-I$ diagrams of the 238 star clusters from \cite{Vansevicius2009}, compared to the grid of star cluster models built in \S\,\ref{sec:grid}.

\begin{figure*}
\centering
\includegraphics[width=180mm]{fig5.pdf}
\caption{\small M31 galaxy star cluster sample of 238 objects from \cite{Vansevicius2009} (circles and stars) and the studied sample of 216 objects (stars) having galactocentric distance over 7 kpc. Panels show: (a) $U-B$ vs $B-V$ and (b) $U-V$ vs $R-I$ diagrams. The model grid (without extinction) used to derive physical parameters of clusters is displayed in the background with indicated reddening vectors following the Galactic extinction law. The continuous line traces PADOVA SSP of $Z=0.008$.}
\label{fig5}
\end{figure*}

Since our grid of models has single metallicity, $Z=0.008$, we attempted to exclude more metal-rich clusters from \cite{Vansevicius2009} sample by only selecting objects with galactocentric distance over 7 kpc. This subsample consists of 216 clusters and is displayed in Fig.\,\ref{fig5}. It is studied with our method described in \S\,\ref{sec:method}, using the Milky Way standard extinction law \citep{Cardelli1989} and distance modulus to M31 of $(m-M)_{0} = 24.47$ derived by \cite{McConnachie2005}.

Figure \ref{fig6} presents the age (panels a, b, c), mass (panels d, e, f) and extinction (g, h, i) of 211 clusters (from the 216 sample) derived by using our method; for five clusters, no model was found within the 3--$\sigma$ around the observed magnitudes. Clusters were studied twice, first with a narrow extinction range $E(B-V)$ = [0.04, 0.5] allowed in the model grid, and second with a wider one, [0.04, 1.0]. The first column (panels a, d, g) compares the age, mass, and extinction derived when a narrow extinction range is allowed vs the \cite{Vansevicius2009} results. The second column (panels b, e, h) compares the age, mass, and extinction derived when a wide extinction range is allowed vs the \cite{Vansevicius2009} results. The last column (panels c, f, i) compares the results obtained with a wide extinction range allowed vs the ones obtained with a narrow extinction range allowed.

\begin{figure*}
\centering
\includegraphics[width=180mm]{fig6.pdf}
\caption{\small M31: age, mass, and extinction of 211 star clusters from \cite[][as V09]{Vansevicius2009}. Panels (a, d, g) show results when a narrow extinction range of $E(B-V)$ = [0.04, 0.5] is allowed in the model grid, compared to V09 values. Panels (b, e, h) show results for a wider allowed extinction range [0.04, 1.0], compared to V09 values. Panels (c, f, i) compare results obtained in a wide extinction range vs the ones obtained in a narrow extinction range. In panels (a, b, c), a density map of Fig.\,\ref{fig3}\,(d) is reproduced to show regions of age-extinction degeneracy. Error bars are computed as described in \S\,2. Six clusters from V09, which have $E(B-V) > 1$, are not shown in panels (g) and (h). Dashed ellipses, numbered ``1'', ``2'', and ``3'', describe majority of associated clusters in the age and extinction panels.}
\label{fig6}
\end{figure*}

In Fig.\,\ref{fig6}\,(a), when clusters are studied in a narrow allowed extinction range, $E(B-V)$ = [0.04, 0.5], the derived ages show the same features as the models studied in Fig.\,\ref{fig3}\,(d), reproduced here in the background of the panel. The degeneracy streaks develop above and below the one-to-one line, perpendicularly to it, and are marked by ellipses numbered ``1'' and ``2'' in Fig.\,\ref{fig6}\,(a). As for the models in the background, the upper degeneracy streak concerns clusters with overestimated age and underestimated extinction, shown by ellipse ``1'' in Fig.\,\ref{fig6}\,(g). In contrast, for the lower degeneracy streak (ellipse ``2''), the age is underestimated and extinction is overestimated. Since the upper streak (ellipse ``1'') is more developed than the lower one (``2''), we interpret that as a result of too narrow an extinction range in the model grid allowed $E(B-V)$ = [0.04, 0.5].

In Fig.\,\ref{fig6}\,(b), when clusters are studied in a wider model extinction range, $E(B-V)$ = [0.04, 1], the upper degeneracy streak is less populated, and the lower one extends, meaning that some of the clusters have overestimated extinction (ellipse ``3''), also shown in panel (h).

If the intrinsic range of extinction for a cluster sample is wide, as in the galaxy M31, then a narrow extinction range of models produces extended an upper streak (ellipse ``1'') and a smaller lower streak (ellipse ``2''), shown in Fig.\,\ref{fig6}\,(a). In panel (b), when the allowed extinction range of models is wide, the upper streak retracts and the lower streak develops. We conclude that in a galaxy with wide extinction range, it is not possible to derive parameters for the clusters affected by age-extinction degeneracy and additional constraints for the extinction are needed; e.g., \cite{Vansevicius2009} used a {\it Spitzer} emission map, to trace the dust lanes of M31 and reduce the age-extinction degeneracy.

Figs. \ref{fig6}\,(d, e) show that for high-mass clusters we obtain lower masses than given by \cite{Vansevicius2009}. Inspecting the metallicity values provided by \cite{Vansevicius2009}, we note that high-mass clusters have metallicities lower than the $Z = 0.008$ used in our model grid. This could be a sign of the metallicity effect on the derivation of physical parameters, and will be investigated in the forthcoming paper.

\subsection{The M33 galaxy star cluster sample}
\cite{SanRoman2009} reports observing of 161 clusters using the {\it HST}, allowing them to evaluate their ages and extinctions by PADOVA \cite{Girardi2002} isochrone fit to their resolved CMDs. Recently, \cite{Ma2012} has used 392 clusters from the compiled catalog of \cite{Sarajedini2007} and images from \cite{Massey2006} to provide $UBVRI$ broad-band integrated photometry data.

There are 40 clusters common to both the \cite{SanRoman2009} and the \cite{Ma2012} catalogs, making them interesting to compare the parameters derived from the resolved method (isochrone fit to CMDs by \cite{SanRoman2009}) and stochastic method, which was applied here to the integrated photometry data of \cite{Ma2012}.

Figure \ref{fig7}\,(a) presents the $U-B$ vs $B-V$ diagram of the 392 clusters from the \cite{Ma2012} catalog compared to the grid of cluster models. The selected 40 star clusters are also indicated. To account for calibration of the $B$ band, as discussed by \cite{Ma2012}, where he compared his photometry to previous works by \cite{Sarajedini2007}, \cite{Park2007}, and \cite{SanRoman2009}, we shifted the $B$ band to brighter magnitudes by 0.1 mag for all clusters. We note that the parameter derivation results do not change when the $B$ band is not taken into account.

\begin{figure*}
\centering
\includegraphics[width=180mm]{fig7.pdf}
\caption{\small M33 galaxy star cluster sample of 161 objects from catalog of \cite{SanRoman2010} (circles and stars) and the studied sample of 40 objects (stars) common to \cite{Ma2012} and \cite{SanRoman2009}. Panels show (a) $U-B$ vs $B-V$, and (b) $U-V$ vs $R-I$ diagrams. The model grid (without extinction) used to derive physical parameters of clusters is displayed in the background with indicated reddening vectors following the LMC extinction law. The continuous line traces PADOVA SSP of $Z=0.008$. The $B$ band magnitudes of clusters were shifted by 0.1 mag.}
\label{fig7}
\end{figure*}

Since the metallicity content of the M33 galaxy is similar to the one of the LMC \citep{Bresolin2010}, we assume that the interstellar extinction law derived by \cite{Gordon2003} for the LMC can be applied to the M33 cluster sample. As in \cite{SanRoman2009}, we adopted the M33 distance modulus of $(m-M)_{0} = 24.54$ \citep{McConnachie2005}.

For the M33, \cite{U2009} derived the extinction of 22 supergiants, which resulted in an extinction distribution centered on $E(B-V)$ = 0.1. \cite{U2009} also used the data of \cite{Rosolowsky2008} to derive $E(B-V)$ values for 58 HII regions, and show in their Fig.\,9 that extinction can be expected to be $E(B-V)$ $\lesssim$ 0.3 for those regions (except for 3 objects), with an average $E(B-V)$ $\sim$ 0.11. A foreground Galactic line-of-sight extinction in the direction of M33 of $E(B-V)$ = 0.04 is estimated from \cite{Schlegel1998} extinction maps.

We studied the M33 clusters using two different allowed extinction ranges for the model grid. The first one is narrow, extending between the foreground Galactic line-of-sight extinction in the direction of M33 up to the higher limit given by \cite{U2009}, i.e, $E(B-V)$ = [0.04, 0.30]. The second one is a wider extinction range, $E(B-V)$ = [0.04, 1.0], such as the one used for the M31 case.

Figure \ref{fig8} presents the age (panels a, b, c), mass (panels d, e, f), and extinction (g, h, i) of the 40 clusters derived by using our method. The results obtained in the narrow extinction range (panels a, d, g) and the wide extinction range (panels b, e, h) are compared to the ones obtained by \cite{SanRoman2009} using isochrone fit to CMDs. Panels (c, f, i) compare the results we obtained with the wide extinction range allowed vs the narrow extinction range allowed. It shows that the difference is relatively small between the two kinds of solution, showing that the age-degeneracy is playing a minor role in this M33 cluster sample.

\begin{figure*}
\centering
\includegraphics[width=180mm]{fig8.pdf}
\caption{\small M33: age, mass and extinction of 40 star clusters common to catalogs of \cite{Ma2012} and \cite[][as S09]{SanRoman2009}. Panels (a, d, g) show results when a narrow extinction range of $E(B-V)$ = [0.04, 0.3] is allowed in the model grid, compared to S09 values. Panels (b, e, h) show results for a wider allowed extinction range of $E(B-V)$ = [0.04, 1.0], compared to S09 values. Panels (c, f, i) compare results obtained in a wide extinction range vs the ones obtained in a narrow extinction range. In panels (a, b, c), a density map of Fig.\,\ref{fig3}\,(d) is reproduced to show regions of age-extinction degeneracy. Error-bars are computed as described in \S\,2.}
\label{fig8}
\end{figure*}

The derived mass does not differ significantly, as seen in Fig.\,\ref{fig8}\,(f). Panels (g, h, i) reveal that only one fourth of the clusters are affected by an increase in extinction when the allowed extinction range is widened. Only some of them have a much higher derived extinction when the wide allowed extinction range is used. These few clusters affected by the age-extinction degeneracy can be seen in panel (c), appearing younger when a wide extinction range is allowed than when the range is narrow, as already seen in M31 cluster sample in Fig.\,\ref{fig6}\,(c).

\section{Conclusions}
We presented a method that aims to derive the age, mass, and extinction of unresolved star clusters by using broad-band photometry and taking the stochastic sampling of stellar masses in clusters into account. We investigated the behavior of the method on a sample of artificial star clusters, in order to trace the different degeneracies between parameters, for different choices of photometric errors and extinction.

These tests allowed us to quantify the age-extinction degeneracy. We demonstrated that for the intrinsically narrow extinction range of the star cluster sample, the age-extinction degeneracy can be resolved even in cases where individual exact extinction values are not known for each cluster.

The age-extinction degeneracy has been observed in the real star cluster sample of \cite{Vansevicius2009} for M31 and to a lesser extent in the one composed of clusters common to the \cite{Ma2012} and \cite{SanRoman2009} catalogs for M33. The physical parameters derived by our method for different extinction ranges in each case, have been compared to the values provided in these studies, showing the impact of the age-extinction degeneracy, especially when the true extinction range of the cluster population is wide.

The M31 star cluster sample from \cite{Vansevicius2009} showed that a true extinction range in this galaxy is wide enough, so that the age-extinction degeneracy is significantly developed, making the parameters difficult to derive, especially for older ages. In such cases, it is preferable to use external constraints on the extinction of individual clusters.

The M33 star cluster sample shows little sign of age-extinction degeneracy, since there is only a small difference between solutions in narrow and wide extinction ranges. The range of allowed extinction of $E(B-V)$ = [0.04, 0.30] gives parameters consistent with the results of isochrone fit to CMDs by \cite{SanRoman2009}.

The follow-up paper will be dedicated to study of the metallicity effects to derive physical parameters of star clusters and to gain a complete understanding of degeneracies introduced by stochastic effects.

\begin{acknowledgements}
We are grateful to the anonymous referee for fruitful comments, which helped to improve the paper. This research was funded by a grant (No. MIP-102/2011) from the Research Council of Lithuania.
\end{acknowledgements}


\section{Introduction}

Dopants in silicon show potential as qubits for solid-state quantum computers~\cite{Kane_Proposal,Skinner_Proposal,Loss_DiVincenzo_Proposal}, with the advantages of scalability as well as the promise of utilizing the
existing semiconductor industry and its processing techniques~\cite{Kane_Proposal,Skinner_Proposal,Loss_DiVincenzo_Proposal,Morton_Nature_Review,Levy_JPCM_Review}.
The theory of Group V dopants
such as phosphorus in silicon~\cite{KohnLuttinger} is useful for describing the quantum nature of the electrons in these systems, as well as for developing schemes to circumvent one of the most challenging aspects of solid-state quantum computers, namely environmental decoherence~\cite{Morton_Nature_Review,Feher_ESR1,Feher_ESR2,Gordon_PRL_58,Witzel_PRB_06,Abe_PRB_10,deSousa_PRB_03,Tyryshkin_PRB_03,Schenkel_APL_06}.
In order to provide a benchmark for such theories and also to use as a starting point for building efficient and accurate tight binding methods, an \textit{ab initio} description of
dopants in silicon is desired.
However, the size of the systems required to describe doped silicon at or near the single-dopant limit is large, making such a description computationally expensive.
In this work, we present large-scale density functional theory calculations for phosphorus-doped silicon supercells with up to 432 atoms.
We make comparisons to other theoretical works~\cite{Sarma_Heitler_London,Carter_Nanotech_11,Kettle_JPCM_04,Wellard_PRB_03}  to determine what can and cannot be captured by approximate or single-electron theories for the doped-silicon systems.

Previous efforts to describe the electronic structure of silicon dopants include effective mass approaches beginning with the work of Kohn and Luttinger~\cite{KohnLuttinger} and continuing with many others~\cite{MacMillen_PRB_84,Kettle_PRB_03,Kettle_JPCM_04,Wellard_JPCM_04,Calderon_PRL_06,Calderon_PRB_07,Hao_PRB_09}, including Fang \textit{et. al.} who perform two-electron Hartree-Fock calculations within effective mass theory~\cite{Fang_PRB_02}.
These efforts include calculating the effects of applied electric and magnetic fields~\cite{Kettle_PRB_03,Calderon_PRL_06} and the coupling of two donors via the exchange interaction~\cite{Kettle_JPCM_04,Wellard_PRB_03,Wellard_JPCM_04,Calderon_PRB_07}.
Tight-binding calculations have also been performed~\cite{Lansbergen_NP_08,Rahman_PRB_09,Rahman_PRB_10}, including a calculation of the quadratic Stark coefficient of the hyperfine interaction which
has reproduced experimentally measured values more accurately than effective mass
theory~\cite{Rahman_PRL_07}.
\textcolor{black}{Two-dimensional layers of dopants in silicon known as $\delta$-layers have been described using density functional theory with compact atomic-orbital basis sets~\cite{Carter_PRB_09,Carter_PRB_09E}, and additional DFT studies~\cite{Carter_Nanotech_11,Carter_PRB_13} evaluated the use of mixed pseudopotentials, which treat the dopant and silicon atoms in the layer using the same core potential, and compared them to all-atom calculations.}
These DFT calculations and a number of additional calculations (see \textcolor{black}{Ref.~\cite{Drumm_NRL_13} and Refs. [30-40] therein}) show a large amount of disagreement for calculated properties such as the valley splitting.
\textcolor{black}{Although some of this disagreement can be attributed to geometrical effects of dopant placement which we will not explore in detail (see instead Refs.~\cite{Drumm_NRL_13}~and~\cite{Carter_Nanotech_11}), the accuracy of the description of dopant electronic structure contributes to these discrepancies.}
Additionally, electrically detected magnetic resonance (EDMR)~\cite{Lo_PRL_11} has called into question a theoretical picture of the scattering of electrons in a two-dimensional electron gas (2DEG) from dopants in silicon~\cite{Ghosh_PRB_92,deSousa_PRB_09}.
An \textit{ab initio} description of a dopant in silicon is therefore useful both as a benchmark and for determining the details of the electronic structure of an isolated dopant which can subsequently be used to calculate more accurate spin-dependent scattering cross sections.
Although these are expensive calculations, we have been able to perform large-scale calculations using the computational resources at Lawrence Livermore National Laboratory.

\section{Methods}
We have used the {\sc Quantum ESPRESSO} suite of programs~\cite{QE-2009} in order to perform density functional theory calculations using a basis set of plane waves.
Face-centered cubic lattices of substitutionally doped silicon were prepared at variable dopant ratios by substituting 1 phosphorus atom in unit cells with 53, 127, 249, and 431 silicon atoms, respectively.  We used the
Perdew-Burke-Ernzerhof (PBE) density functional~\cite{PBE1,PBE2} with an ultrasoft pseudopotential for phosphorus (P.pbe-n-van.UPF from Ref.~\cite{qewebsite}) and a norm-conserving
pseudopotential that was calculated using FHI98PP\cite{Fuchs99} for silicon.
A plane-wave energy cutoff of 65~Ry. was chosen based on convergence of the total energy and pressure in the smallest supercell.
K-space sampling was performed using a Monkhorst-Pack~\cite{MonkhorstPack} grid of 8$\times$8$\times$8 (54 atom), 6$\times$6$\times$6 (128 atom), 4$\times$4$\times$4
(250 atom), and 2$\times$2$\times$2 (432 atom) grid points.
As a reference for the energy and properties of bulk silicon, a two-atom silicon cell was used with a grid of 20$\times$20$\times$20 k-points.
In this cell, one silicon is placed at the origin and another is placed at the point $(\frac{a}{4},\frac{a}{4},\frac{a}{4})$, where $a$ is the lattice constant.
The basis vectors of the FCC cell are $(\frac{a}{2},\frac{a}{2},0)$, $(\frac{a}{2},0,\frac{a}{2})$, and $(0,\frac{a}{2},\frac{a}{2})$; the two silicon atoms are repeated at every integer multiple of the basis vectors.
The lattice constant of the doped supercells was determined as multiples of the 5.46 \AA~lattice constant for bulk silicon computed with the pseudopotential used in this study.
\textcolor{black}{
The phosphorus dopants repeat along the directions of the basis vectors.
The directions of the plots given below were chosen to be orthogonal to each other.
Therefore, for the (001) plane, first the [100] and [010] basis vectors were chosen, and then the Gram-Schmidt procedure was used to find the orthogonal direction ([-1 2 0]).
The angular bracket notation is used in the plots because the choice of direction is somewhat arbitrary due to the FCC symmetry of the system.
}
Geometric minimization of the total energy was performed for all doped systems studied here. In the calculations with $N<432$ atoms, all of the atoms were
allowed to move during the simulation, while for $N=432$, the atoms at the edges of the supercell were frozen in order to make better comparisons to the bulk silicon system.  The
volume of each supercell was constant during the geometric optimization.
These computations required 600 processors on the Ansel supercomputer at Lawrence Livermore National Laboratory, a 324 node system using the Intel Xeon architecture rated at a peak of 43.5 TFLOP/s~\cite{llnlcomputing}.
A single point energy calculation for the 432-atom cell took about 19 CPU hours to complete.

In Section~\ref{potentialsection}, we use the doping potentials as a measure of the electronic environment of the dopants.
The electric potential for each cell is obtained by adding the nuclear term ($V_{nuc}$), coulombic Hartree term ($V_{Hartree}$), and PBE exchange-correlation term ($V_{xc}$) at the converged electronic density $\rho _{conv}$,
\begin{equation}
\label{cellpotential}
V = V_{nuc} + V_{Hartree} [\rho _{conv}] + V_{xc} [\rho _{conv}].
\end{equation}
The density $\rho _{conv}$ is that which minimizes the energy of the system according to Eq.~(\ref{cellpotential}).
The doping potential is then obtained by subtracting the potential obtained for the undoped cell from that of the doped cell,
\begin{equation}
\label{dopingpot}
V_{doping} = V_{doped} - V_{undoped}.
\end{equation}
Similarly, the doping density discussed in Section~\ref{densitysection} is obtained by subtracting $\rho _{conv}$ of the undoped cell from that of the doped cell.
Finally, the spin densities discussed in Section~\ref{spinsection} are obtained by subtracting the spin-up and spin-down portions of the density $\rho _{conv}$ of a single cell.

In order to estimate the exchange coupling (J), we use the DFT broken spin symmetry states~\cite{Rudra_JCP_06,Ruiz_JCC_11,Yamaguchi_CPL_79,Noodleman_JACS_85,Noodleman_JCP_81},
\begin{equation}
\label{jcoupling}
\frac{J}{2}=E_{BS}-E_{HS},
\end{equation}
where $BS$ and $HS$ denote, respectively, the broken symmetry ground state and high-spin state in which the spin density is constrained to give a spin of $S=\frac{1}{2}$ in each unit cell.
This procedure empirically has given exchange couplings with a satisfactory degree of accuracy for molecular systems~\cite{Rudra_JCP_06,Ruiz_JCC_11}, and the exchange couplings it calculates can be thought of as the spin-coupling parameter $J$ in both Ising and Heisenberg models of the spin interactions in the periodic system~\cite{Illas_TCA_00,Dai_JCP_03}.

\section{Results and Discussion}
The Bohr radius of phosphorus in silicon is estimated to be about 2.5 nm~\cite{KohnLuttinger}.
In order to reach the single dopant limit, the phosphorus atoms in a silicon matrix would have to be spaced a minimum of $\sim5$~nm apart.
This would require a supercell of tens of thousands of atoms and is outside the realm of feasibility for most DFT calculations.
Here, we have therefore assessed the energies and properties of doped silicon as a function of the doping fraction approaching this limit, in order to develop a better
understanding of how explicitly treating the phosphorus donor atom and lattice relaxation effects the electronic structure.
\textcolor{black}{In Table~\ref{energytable}, we
show the convergence (as a function of the size of the supercell) of
the geometrically relaxed and unrelaxed formation energy,
\begin{equation}
\label{formenergy}
E_{doped}-(E_P+(N_{atom}-1)E_{Si}) - (E_{undoped}+N_{atom}E_{Si}),
\end{equation}
where $E_{doped}$ is the energy of the doped supercell at either the geometrically optimized (relaxed) or bulk silicon (unrelaxed) geometry, $E_P$ is the energy of an isolated phosphorus atom, $E_{Si}$ is the energy per atom of bulk silicon, $E_{undoped}$ is the energy of the undoped supercell of $N_{atom}$ atoms, and the last term in parentheses in Eq.~(\ref{formenergy}) is included to account for using approximate energy cutoffs and k-point grids for the larger supercells.
\begin{table}[htp!]
\color{black}
\caption{Convergence of the energy as a function of unit cell size.
For each unit cell from 54 to 432 atoms, we give the unrelaxed formation energy (column 2), the relaxation energy (column 3), and the relaxed formation energy (column 4)
\label{energytable}}
\begin{ruledtabular}
\begin{tabular}{cddd}
\multicolumn{1}{c}{Unit cell}&
\multicolumn{1}{c}{Formation Energy (unrelaxed)}&
\multicolumn{1}{c}{Relaxation Energy}&
\multicolumn{1}{c}{Formation Energy (relaxed)}\\
&
\multicolumn{1}{c}{(eV)}&
\multicolumn{1}{c}{(meV)}&
\multicolumn{1}{c}{(eV)}\\
\hline
Si$_{53}$P&-3.1519&-30.0512&-3.1820\\
Si$_{127}$P&-3.1816&-27.8154&-3.2095\\
Si$_{249}$P&-3.1872&-32.1070&-3.2193\\
Si$_{431}$P&-3.2322&-32.3320&-3.2646
\end{tabular}
\end{ruledtabular}
\end{table}
The geometric relaxation energy of the doped system (column 3 of Table~\ref{energytable}) is similar (around 30~meV) for each of the cells.
The formation energies increase in magnitude as a function of cell size, suggesting that the presence of the defect is reducing the strain in the lattice.
The relaxed formation energy increases in magnitude by 82.6 meV from the 54 atom cell to the 432 cell and by 45.3 meV from the 250 to the 432 atom cell.
This suggests that the effects of lattice relaxation and changes in the electronic structure will be important for the donor electron dynamics in phosphorus doped silicon where
the isotropic doping fraction is $\sim 0.2$~\% or higher.}  In systems with a lower doping fraction, these effects may still play a role, but we cannot make a definitive statement on this issue
since the changes in energy exhibit a strong dependence on the system size up to the largest system studied here.  In another recent study\textcolor{black}{~\cite{Drumm_NRL_13}}, the energy gained by relaxing a cell with a monolayer of phosphorus
dopants was found to be of a similar magnitude.  Whether the lattice relaxation is important to EDMR readout schemes depends on how the relaxation affects the scattering dynamics of electrons, which is related to the doping potential.

In the next section, we make comparisons of the doping potentials
for a donor phosphorus atom in silicon, both with and without the effects of geometrical relaxation of the lattice.

\subsection{Doping potential\label{potentialsection}}
We define the doping potential as the doped cell potential minus the bulk silicon potential (Eqs.~(\ref{cellpotential})~and~(\ref{dopingpot})).
The doping potential shows how the electronic environment surrounding the dopant differs from that of bulk silicon.
These potentials also provide input for calculations of the cross sections of electron scattering at the dopants~\cite{Lordi_Scattering_PRB_10,SCARLET},
which are largely determined by integrals of the doping density~\cite{Lordi_Scattering_PRB_10}.
By calculating the scattering of conduction electrons confined in a two-dimensional layer located at a given distance from the (001) plane, a connection can be made with electrically detected magnetic resonance (EDMR) schemes~\cite{Ghosh_PRB_92,Xiao_Nature_04,Lo_PRL_11,vanBeveren_APL_08,Lo_APL_07} used to measure the dopant spin state.
Doping potentials calculated using atomistic DFT can also be used to parameterize new tight binding models, or effective single-electron models which more accurately reproduce the effects of the dopant electronic structure than standard effective mass models.

The doping potentials  are shown for the 54 and 432 atom cells in Figs.~\ref{potentialfigure}~and~\ref{relaxedpotfigure}.
\textcolor{black}{In Fig.~\ref{potentialfigure}, the potentials are given for both the doped and undoped cells at geometrically unrelaxed (Figs.~\ref{m54u}~and~\ref{m432u}) and relaxed (Figs.~\ref{m54}~and~\ref{m432}) geometries.  In Fig.~\ref{relaxedpotfigure}, the effects of geometric relaxation on the doping potentials is given.}
\begin{figure}[htp!]
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI53P_Doping_Bulk_Cax2}\label{m54u}}
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI431P_Doping_Bulk_Cax2}\label{m432u}}\\
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI53P_Doping_Optimized_Cax2}\label{m54}}
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI431P_Doping_Optimized_Cax2}\label{m432}}
\caption{(Color online) The potential difference (eV) in the (001) plane between doped and undoped silicon is shown for cell sizes of 54~(\ref{m54u}~and~\ref{m54}) and 432 atoms~(\ref{m432u}~and~\ref{m432}).
An area of 93 $\times$ 55 \AA  is shown for in all figures, in which the 54 and 432 atom cells repeat about 44 and 11 times, respectively.
For the 54 and 432 atom cells, parts of 51 and 14 dopants, respectively, can be seen in the figures.
Doping potentials are shown for cells at the bulk geometry~(\ref{m54u}~and~\ref{m432u}) and optimized geometries~(\ref{m54}~and~\ref{m432}).
The contour lines are drawn every 1 eV except between -2 and 2 eV, where they are drawn every 100 meV.
The color axis in these figures is between -1 and 1 eV in order to highlight the effects at this energy scale, while the range of doping energies is between -9 and 12 eV, with the larger values occuring near the dopant nucleus (see Fig.~\ref{potentialcufigure}).
\label{potentialfigure}}
\end{figure}
\begin{figure}[htp!]
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI431P_RelaxingDoping_Cax2}\label{m432rd}}
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI431P_RelaxingDoping_Cax3_Z}\label{m432rdz}}
\caption{(Color online) The potential difference (eV) in the (001) plane between doped silicon in its optimized geometry and bulk silicon for the 432 atom cell in the bulk geometry.
We refer to this as the ``relaxing and doping potential.''
In Fig.~\ref{m432rd}, 14 phosphorus dopants are visible, while Fig.~\ref{m432rdz} shows a close-up of the region around one dopant.
The contours and color axis in Fig.~\ref{m432rd} are the same as in Fig.~\ref{potentialfigure}.
The close up in Fig.~\ref{m432rdz} uses the same contours, but the color axis is between -2 and 2 eV.
\label{relaxedpotfigure}}
\end{figure}
In the (001) plane, the doping potentials for the 54 atom cell (left panels) can be seen to overlap.
Thus, while there are areas in which the doping potential goes to zero, indicating a return to bulk silicon behavior, the dopants are largely connected by regions of potential greater than 100 meV.
In the 432 atom cell at the bulk geometry, the doping potentials nearly go to zero between the dopants.
However, there are small areas of non-zero potential still connecting dopants, suggesting that the single dopant limit has still not been reached.
These potential connections between dopants are slightly exaggerated at the relaxed geometry.
In both cells, however, the potential region directly around the dopants is similar, with a region exceeding 8 eV directly around the dopant and three 1 eV lobes about 120 degrees apart from each other.
\textcolor{black}{These lobes are evident even as geometric relaxation effects are included (Fig.~\ref{relaxedpotfigure}).  Additional oscillations are evident in the potential when geometric relaxation effects are included as a result of the position of the nearby silicon nuclei, although the pattern of oscillations near the dopant are similar with and without relaxation effects.}
In addition to the lobes near the nucleus, there are additional lobes at an energy scale of -700 meV which extend away from the nuclei in space.
The trigonal symmetry results from taking a two-dimensional cut through the three dimensional structure onto the (001) plane.
When viewed in three dimensions, these lobes can be seen to arise from the tetrahedral nature of the bonding in the silicon cell.
Three-dimensional potential isosurfaces at -700, -600, and -150 meV are shown in Fig.~\ref{3dpotentialfigure}.
\begin{figure}[htp!]
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{si432p_n700_redder}\label{3d700}}
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{si432p_n600_redder}\label{3d600}}\\
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{si432p_n150}\label{3d150}}
\caption{(Color online) Doping potential isosurfaces are shown for -700~(\ref{3d700}), -600~(\ref{3d600}), and -150~(\ref{3d150}) meV.
The isosurfaces are shown in slightly transparent blue, and the phosphorus donors are shown in red.
The silicon atoms are not pictured in Figs.~\ref{3d700}~and~\ref{3d600}, and they are shown for reference in Fig.~\ref{3d150}.
\label{3dpotentialfigure}}
\end{figure}
\textcolor{black}{
To address the question of whether the tetrahedral doping potential configuration arises from the FCC structure of the periodic cell images, we have also performed calculations for rectangular cells with an aspect ratio of two to one.
The dopant potentials for these cells (not pictured) exhibit tetrahedral doping potentials similar to those seen in Figs.~\ref{potentialfigure},~\ref{relaxedpotfigure},~and~\ref{3dpotentialfigure}.
This suggests that the salient features of the doping potential such as the tetragonal lobes are not dependent on the geometry of the cell arrays.
}
The non-spherical nature of the doping potential will be important for calculation of electron-dopant scattering cross sections.
Effective mass models of doped silicon assume a spherical Coulomb potential, while some tight binding models use a spherical, Coulomb-like doping potential~\cite{NEMO3D,Rahman_PRL_07,Rahman_PRB_10,Rahman_PRB_11} and allow the surrounding atoms' electronic structure to adjust according to this potential.
The anisotropic doping potentials and cell geometries calculated here could  be used as alternative parameterizations for tight binding models, and they can also be used to calibrate the resulting potentials and densities calculated by tight binding models.

Fig.~\ref{potentialcufigure} shows a closer view of the doping potential in the region of the dopant for the 432 atom cell with optimized geometry.
\begin{figure}[htp!]
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_Doping_Optimized_Cax1_Z}}
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_Doping_Optimized_Cax1_ZZ}}
\caption{(Color online) The potential difference (eV) between doped and undoped silicon is shown for the 432 atom cell at the optimized geometry of the doped cell for the region close to a  dopant.
\textcolor{black}{Different zooms are shown in panels (a) and (b).}
The doping potential can be seen to oscillate in this region as opposed to showing s-orbital character as predicted by effective mass theory.
Additionally, there are a number of features which are not spherically symmetric.
The contour lines are the same as in Fig.~\ref{potentialfigure}, while the color axis has been expanded to the range of -9 to 12 eV.\label{potentialcufigure}}
\end{figure}
In the effective mass picture~\cite{KohnLuttinger}, the dopant wavefunction has s-orbital character.
In Fig.~\ref{potentialcufigure}, the potential in the region of the dopant can be seen to oscillate as a function of the orientation.
The oscillations are due to interactions with electrons in the shells below the valence shell.
If these calculations were performed without using pseudopotentials, which reduce oscillations from core electrons and replace them with a smooth potential, the potential would most likely oscillate to an even greater degree.
These oscillations, as well as those visible in the optimized geometries of Fig.~\ref{potentialfigure} due to the silicon lattice distortions, represent qualitative differences between DFT and effective mass theory.

In Ref.~\cite{Carter_Nanotech_11}, Carter \textit{et. al.} use mixed pseudopotentials in order to estimate the potential as a function of the distance from a layer of dopants.
In Fig.~\ref{dopinglayerfigure}, we plot the doping potential (doped minus bulk) as a function of distance from the (001) layer of dopants.
\begin{figure}
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_DopingD2_Optimized_Cax2}\label{dpsi}}\\
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_DopingD2_1D}\label{dpmr0}}
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_RelaxingDopingD2_1D}\label{rdpmr0}}
\caption{(Color online) The doping potentials ($V_{Si_{431}P}-V_{Si}$) as a function of distance from the doping layer.
Fig.~\ref{dpsi} gives a contour plot with the x-direction representing an intralayer direction and the y-direction representing an interlayer direction.
The contours and color axis are the same as in Fig.~\ref{potentialfigure}.
Fig.~\ref{dpmr0} shows the potential as a function of distance from the doping layer along the line between dopants.
The inset shows an enlargement for the energy range -200 to 200 meV, which is most relevant to scattering by a two-dimensional electron gas.
Fig.~\ref{rdpmr0} shows the potential with geometric relaxation included.\label{dopinglayerfigure}}
\end{figure}
In Fig.~\ref{dpsi}, a full two-dimensional cut through the potential is given in the plane perpindicular to the (001) plane, and in Figs.~\ref{dpmr0}~and~\ref{rdpmr0} the potential is shown in a slice of this plane which connects two dopants.
In contrast to Ref.~\cite{Carter_Nanotech_11}, where no structure is evident, marked structure is seen in Fig.~\ref{dopinglayerfigure}.
In Figs.~4~and~7 of Ref.~\cite{Carter_Nanotech_11}, the mixed pseudopotential doping potentials are much smoother than in Fig.~\ref{dpmr0}, especially in the region around the dopant.
In Fig.~\ref{dpmr0}, there is a significant amount of structure in the potential near the dopant itself.
Minor effects of the silicon atoms in the next layer of the crystal are also evident in Fig.~\ref{rdpmr0} when the effects of geometric relaxation are included.
It is important to note that these results are for doping densities near the single dopant limit: a study of the effect of a $\delta$-layer of dopants would require very large cells which would likely have thousands of atoms.
\textcolor{black}{Additionally, the dopant potentials in Ref.~\cite{Carter_Nanotech_11} are plane-averaged, while we have plotted straight point potentials.
However, averaging does not eliminate the structure in our potentials, but instead reduces the peak potential relative to the somewhat noisy structure of the atomic lattice.}

\subsection{Constrained spin calculations\label{spinsection}}
In the effective mass model~\cite{KohnLuttinger,MacMillen_PRB_84,Kettle_PRB_03,Kettle_JPCM_04,Wellard_JPCM_04,Calderon_PRL_06,Calderon_PRB_07,Hao_PRB_09},
the dopants are treated as effective hydrogenic one-electron systems, with spin due to the additional dopant electron.
Density functional theory gives a more detailed treatment of the many-body problem,

allowing dopant electrons to couple to core electrons on the same dopant, to electrons in the silicon atoms, and
most importantly, also allowing dopant electrons to couple with each other at low densities.
In order to compare the DFT results with the frequently used single-electron picture~\cite{Morton_Nature_Review,Feher_ESR1,Feher_ESR2,Gordon_PRL_58,Witzel_PRB_06,Abe_PRB_10,Tyryshkin_PRB_03,Schenkel_APL_06}, we calculated the spin densities (spin up density minus spin down density) of each doped cell.

Fig.~\ref{spindensityfigure} shows spin densities for DFT calculations in which the total spin of the individual unit cells was constrained to be $S=\frac{1}{2}$ a.u. (one bohr magneton) \textcolor{black}{along the z-axis}.
\begin{figure}[htp!]
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI53P_SpinDensity_Cax1}\label{s54}}
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI127P_SpinDensity_Cax1}\label{s128}}\\
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI249P_SpinDensity_Cax1}\label{s250}}
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI431P_SpinDensity_Cax1}\label{s432}}\\
\subfigure[]{\includegraphics[height=2in,width=2in,angle=0.0]{SI431P_SpinDensity_Cax1_Z}\label{s432cu}}\\
\caption{(Color online) DFT calculations of spin densities ($\rho _\uparrow -\rho _\downarrow$) of the 54, 128, 250, and 432-atom unit cells in a 15 $\times$ 15 \AA~box centered on a dopant, with
the spin constrained to $S=\frac{1}{2}$ \textcolor{black}{along the z-axis} in each unit cell.
A 2.5 $\times$ 2.5 \AA~close up of the dopant area of the spin density for the 432-atom cell is also shown~(\ref{s432cu}).\label{spindensityfigure}}
\end{figure}
\textcolor{black}{
The spin density is $\rho _\uparrow -\rho _\downarrow$, where $\rho _\uparrow$ and $\rho _\downarrow$ represent the density of electrons with spin aligned or anti-aligned with the z-axis, respectively.
In all of the previous sections, the calculations were unconstrained in the spin degree of freedom and \textcolor{black}{since the ground state of silicon is not magnetic, these all converged to an average spin of $S=0$, meaning the unpaired electron points in random directions with none preferred over any others.}
While the electron donated by the phosphorus atom does indeed have a spin of $\frac{1}{2}$, in the ground state it can be considered as being in a superposition of its aligned and anti-aligned states, leading to a total spin of $S=0$.
In the following spin-constrained calculations, an effective external field is added in the form of an energy penalty,
\begin{equation}
\label{energypenalty}
E_{total} = E_{DFT}+\lambda \left( \rho _\uparrow - \rho _\downarrow - \frac{1}{2}\hat{z} \right).
\end{equation}
}
In addition to giving
an indication of spin ordering due to the coupling between dopant atoms in the system with respect to cell size, the local spin density also gives a qualitative picture of where the ``additional'' electron provided by the phosphorus dopant is located
when the total magnetization of the system is locally constrained.
In the 54 and 128 atom unit cells, the spin density clearly shows interactions between the dopants in the form of areas of high density between dopants.
The 250 and 432 atom cells show comparatively less localization of spin density between the dopant atoms,
but there are still areas of enhanced spin density connecting the dopants along lines in the (001) plane.
These areas are fading for the 432 atom cell but not completely absent.
Note also that the maximum spin density is decreasing from the 54 to the 250 atom cell, but remains approximately constant between the 250 and 432 atom cells.
These results give an effective one-electron picture of how modulations in the spin density are affected by the distance between dopant atoms and provide insight on
the behavior of the effective one-electron wavefunction in this system.

An exchange coupling between donors can be estimated using density functional theory according to Eq.~(\ref{jcoupling})\cite{Rudra_JCP_06,Ruiz_JCC_11,Yamaguchi_CPL_79,Noodleman_JACS_85,Noodleman_JCP_81,Illas_TCA_00,Dai_JCP_03}.
The quantity $J$ provides an estimation of the spin couplings between donor qubits which may be used to apply two-qubit gates.
The exchange couplings calculated for the different size supercells are given in Table~\ref{exchangetable}.
\begin{table}[htp!]
\caption{The exchange parameter as a function of cell size.\label{exchangetable}}
\begin{ruledtabular}
\begin{tabular}{cdd}
\multicolumn{1}{c}{Unit Cell}&
\multicolumn{1}{c}{Distance between dopants (\AA)}&
\multicolumn{1}{c}{Exchange coupling (meV)}\\
\hline
Si$_{53}$P&11.62&-125.1\\
Si$_{127}$P&15.49&-83.6\\
Si$_{249}$P&19.36&-64.4\\
Si$_{431}$P&23.23&-49.8
\end{tabular}
\end{ruledtabular}
\end{table}
As the dopant density decreases, we see that the magnitude of the exchange coupling also decreases.
In Refs.~\cite{Kettle_JPCM_04}~and~\cite{Wellard_JPCM_04}, the exchange coupling of a two-donor system is calculated for dopants at much greater distances than in this work, using a Heitler-London approximation with variable alignment of the longitudinal and transverse Bohr radius of the dopant relative to the inter-dopant direction, respectively.
These works, together with Refs.~\cite{Calderon_PRB_07}~and~\cite{Koiller_PRL_01}, have concluded that oscillations in the exchange coupling would make it difficult to control a quantum information system which attempts to exploit this coupling.
Although donor spacings studied here are small compared to the previous studies of a pair of phosphorus donor atoms in bulk Si in Refs.~\cite{Kettle_JPCM_04}~and~\cite{Wellard_PRB_03}, we may nonetheless make a comparison of our results with
these in Fig.~\ref{exchangefigure}.
For the three largest systems studied here, the behavior of the exchange coupling with respect to donor spacing ($r$) is fit well by a single exponential decay, $J(r)=-323.0\exp(-r/12.0)$.
The exchange coupling at the distances studied using the
Heitler-London approximation is systematically larger than the corresponding value extrapolated from our fitted decay.  The source of the decrease in magnitude of the DFT estimated exchange coupling is likely due to oscillations in the donor electron densities
which are present in the DFT calculation, but which are not included in the models used in
Refs.~\cite{Kettle_JPCM_04}~and~\cite{Wellard_PRB_03}.  We note that due to the relatively small spacing between donors and the isotropic distribution of donors in this work, the correlations between donor
atoms are stronger than what is modeled in the studies of isolated pairs of atoms separated by a large distance.  However, if in actual devices the donor atoms are not evenly distributed and well separated, the magnitude of the exchange coupling may also be
decreased
with respect to what is predicted by the Heitler-London model.   Extending our studies to larger systems with anisotropic doping will allow for a more direct comparison with studies based of isolated donor pairs.
\begin{figure}
\begin{center}
\includegraphics[height=3in,width=3in,angle=0.0]{ExchangePlot}
\end{center}
\caption{Magnitude of the zero-field exchange coupling calculated by DFT in this work (circles), and in the Heitler-London approximation with effective mass theory wavefunctions in Ref.~\cite{Kettle_JPCM_04}~(triangles., Bohr radius of 2.381 nm), and in Ref.~\cite{Wellard_PRB_03}~(squares, Bohr radius of 1.368 nm).  Data for the three largest systems studied here are fit well by a single exponential decay in the spacing between
donors (solid line.)  The magnitude of the exchange coupling predicted by the fitted decay function for isotropic doping is much lower than that calculated in Refs.~\cite{Kettle_JPCM_04}~or~\cite{Wellard_PRB_03} using the Heitler-London approximation for
an isolated pair of P atoms in bulk Si (inset.)\label{exchangefigure}}
\end{figure}

\subsection{Doping density\label{densitysection}}
In Ref.~\cite{Sarma_Heitler_London}, wavefunctions of phosphorus dopants were calculated using effective mass theory, in which the Bloch functions of silicon where taken directly from a density functional theory calculations.
Doping densities (density of doped cell minus undoped cell) were presented in Fig. 2 of Ref~\onlinecite{Sarma_Heitler_London}.
In Fig.~\ref{dopantdensityfigure}, we present the doping density,
\begin{equation}
\rho _{doping} = \rho _{doped}-\rho _{undoped},
\end{equation}
calculated for the 432 atom cell using a full atomistic DFT treatment, where each density $\rho$ is calculated from the Kohn-Sham orbitals $\phi _i$ as
\begin{equation}
\rho = \sum _i \vert \phi _i \vert ^2.
\end{equation}
\begin{figure}[htp!]
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_DopingDensity_Optimized_Cax2}}
\subfigure[]{\includegraphics[height=3in,width=3in,angle=0.0]{SI431P_DopingDensity_Optimized_Cax2_Z}\label{ddzoom}}
\caption{(Color online) The doping density ($\rho _{Si_{431}P}-\rho _{Si}$) in the (001) plane.
The region of the density near the dopant is zoomed in panel~\ref{ddzoom}.
The contours are drawn every 0.01 units, except between -0.1 and 0.0 where they are drawn every 0.002 units.
The color axis is set between -6 and 7 $\times$ 10$^{-3}$ units, while the density varies between -6 $\times$ 10$^{-3}$ and 0.35 units with the larger values near the nucleus.\label{dopantdensityfigure}}
\end{figure}
The current DFT doping densities show a different distribution in the (001) plane than those presented in Ref.~\cite{Sarma_Heitler_London}.
In particular, the present densities are more circular around the dopant in this plane.
The oscillations which are not included in the effective mass calculations are also evident in the immediate vicinity of the dopant.
These oscillations are about two orders of magnitude less than the doping density in the vicinity of the dopants.
Lobes similar to those seen in the doping potentials in Section~\ref{potentialsection} are also apparent in the doping densities.
Finally, we note that the doping density shows evidence of the interaction between dopants, resulting in finite electron density between dopants.

\section{Conclusion}
We have presented density functional theory calculations for
silicon doped with a single phosphorus atom, for systems with up to 432 atoms in a cell.
A detailed knowledge of the electronic structure of doped silicon is necessary for the implementation of spin-based qubits in silicon~\cite{Kane_Proposal,Skinner_Proposal,Loss_DiVincenzo_Proposal,Levy_JPCM_Review}.
We have calculated non-spherical electron densities and doping potentials which allow for a level of microscopic description beyond effective mass and tight binding theory.
In comparison to previous calculations of dopant electronic structure~\cite{KohnLuttinger,MacMillen_PRB_84,Kettle_PRB_03,Kettle_JPCM_04,Wellard_JPCM_04,Calderon_PRL_06,Calderon_PRB_07,Hao_PRB_09,Fang_PRB_02,Lansbergen_NP_08,Rahman_PRB_09,Rahman_PRB_10,Carter_PRB_09,Carter_PRB_09E,Carter_Nanotech_11}, we have found an unprecedented level of structure in the doping potentials (Fig.~\ref{potentialfigure}) and densities/wavefunctions (Fig.~\ref{dopantdensityfigure}).
\textcolor{black}{Due to the oscillatory nature of doping potentials, the exchange coupling between qubits obtained by extrapolating our results to smaller distances was found to be less than estimates based on the Heitler-London approximation, although further calculations at larger dopant separations are required to confirm this result.}

These calculations have many potential uses.
Such detailed microscopic calculations will allow more accurate and detailed device simulations than are currently possible.
By understanding the effects of modulations in the doping density including effects of both the spin density as well as the doping potential, allows now calculations which probe the readout properties of the qubits, especially using techniques such as EDMR~\cite{Ghosh_PRB_92,Xiao_Nature_04,Lo_PRL_11,vanBeveren_APL_08,Lo_APL_07}.
Additionally, alternative qubits such as excited-state dopants~\cite{Stoneham_JPCM_03} or charged dopant qubits~\cite{Sellier_PRL_06,Calderon_PRB_10,Rahman_PRB_11,Taniguchi_SSC_76} can now be explored with this accurate picture of the electronic structure.
The doping potentials calculated here can provide the starting point for effective one-electron calculations of a dopant electron wavefunction, possibly deformed by some electrostatic gate potential.
They also can guide design of multiple qubit devices by providing effective Hamiltonians or potentials for multiple-qubit geometries.
Finally, the doping potentials provide input for scattering calculations, including calculations in which the current-carrying electrons are confined to a two-dimensional plane to model electrical readout schemes for silicon quantum computation.
The spin densities we have calculated here can also be used to compare with a single-electron picture and to determine the density of the electron which is donated by the phosphorus.

In the future, we plan to look at systems which more accurately represent the experimental devices.
This will require looking at the effect of the silicon dioxide interface and defects at this interface, a considerably more computationally intensive task.
They may also be extended to the calculation of parameters related to the hyperfine interaction.
These DFT calculations may also be coupled with calculations of the spin-dependent scattering~\cite{SCARLET,Lordi_Scattering_PRB_10} of the two-dimensional electron gas as well as quantum control calculations for the implementation of quantum logic operations~\cite{Viola_PRL_99,KBW_Noise_PRA_06,KBW_Noise_PRA_08,KBW_QND_PRB_08,KBW_arXiv_10}.

\begin{acknowledgments}
The authors wish to acknowledge insight due to conversations with Dr. Vincenzo Lordi (Lawrence Livermore National Laboratory).
Additionally, we wish to thank our experimental collaborators in the group of Dr. Thomas Schenkel (Lawrence Berkeley National Laboratory), especially Dr. Cheuk Chi Lo and Christoph Weis, for helpful conversations and direction.
This work was supported by the UC Lab Fees Research Program under a grant to the University of California, Berkeley and Lawrence Livermore National Laboratory.  It was also performed in part (HDW) under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.
\end{acknowledgments}

\bibliographystyle{apsrevM}