Mapping Stellar Surfaces. I. Degeneracies in the Rotational Light-curve Problem

In general, inferring all of the properties of a stellar surface from its light curve alone is not only difficult but formally impossible. To understand why, consider an expansion of the stellar surface intensity in the spherical harmonic basis out to arbitrary order ¹² (see Figure 2). Assuming (for the moment) that the star rotates about an axis that points up along the page, the observed light curve may be expressed as a weighted sum of the disk-integrated intensity of each of the spherical harmonics as they rotate about that same axis (Luger et al. 2019). However, not all spherical harmonics will contribute to the full light curve, as many (in fact, most) of the spherical harmonics are perfectly antisymmetric about the equator. This is the case for the l = 1, m = −1 harmonic, which integrates to zero regardless of the phase at which it is viewed. The same is true, in fact, for all other harmonics of order m = −1 and (perhaps less obviously) all harmonics with odd l = 3, 5, 7, .... Furthermore, there exist many linear combinations of spherical harmonics that similarly integrate to zero at all rotational phases. Together, these modes constitute the null space of the problem: the set of modes on the surface that do not contribute to the observed light curve and therefore cannot be probed from photometry.

For rotational light curves, the vast majority of the surface modes lie in the null space. To show this, we will make use of the fact that we can express the vector of K corresponding to observed fluxes f (i.e., the light curve) as a linear operation on the vector of N spherical harmonic coefficients y (Luger et al. 2019):

$$\begin{eqnarray}&&{\boldsymbol{f}}={\bf{1}}+{\boldsymbol{ \mathcal A }}\,{\boldsymbol{y}},\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{ \mathcal A }}$ is the (K × N) design matrix of the transformation, whose columns describe how each of the N components in the spherical harmonic basis contributes to each of the K points in the light curve. ¹³ Even though we are explicitly choosing the spherical harmonics as the basis in which we describe the stellar surface, Equation (1) is quite general and applies to any basis that is linearly related to the flux. For instance, y could instead describe the intensity deficits in each of the N pixels of a gridded stellar surface, in which case ${\boldsymbol{ \mathcal A }}$ would be the matrix of pixel visibilities that describe how to sum each of the pixels to obtain the observed light curve.

The size of the null space is called the nullity, and it is equal to N − R, where N is once again the number of coefficients describing the stellar surface and R is the rank of the flux operator ${\boldsymbol{ \mathcal A }}$. The rank R is the number of linearly independent columns in ${\boldsymbol{ \mathcal A }}$, which can be computed numerically using any standard linear algebra package. It is equal to the number of independent components that can be measured given an observation of f .

Figure 3 shows the rank and nullity of the flux operator as a function of the resolution of the surface map (quantified as the spherical harmonic degree l of the expansion). The orange curve is the number of spherical harmonic coefficients needed to represent a surface map up to degree l and is equal to $N={\left(l+1\right)}^{2}$. The blue line shows the rank of the flux operator ${\boldsymbol{ \mathcal A }}$ as a function of l, which is equal to

$$\begin{eqnarray}\begin{array}{rcl}R & = & 4\left\lfloor \displaystyle \frac{l}{2}\right\rfloor +1\\ & \approx & 2l+1.\end{array}\end{eqnarray} \tag{ 2 }$$

The nullity is simply the difference between N and R.

Zoom In Zoom Out Reset image size

The most striking feature in Figure 3 is how quickly the two curves diverge as l increases. What this means for the mapping problem is that the number of surface modes—the total information needed to represent a surface map at some resolution—grows much more quickly than the number of independent degrees of freedom in the light curve. At all but the lowest resolutions, there are always more features in the null space than components one can measure in the light curve, a difference that grows quadratically with l. This means that although a light curve can tell us some information about a stellar surface on very large scales, the amount of information it tells us quickly decreases for smaller scales and all but vanishes for the smallest surface features.

It is also worth noting the piecewise nature of the rank as a function of l; increasing the degree of the expansion from even l to odd l does not increase the rank of the flux operator. Put another way, odd spherical harmonic modes with l > 1 are always in the null space of the light curve. We will return to this point below.

In the Appendix, we show that it is straightforward to construct linear operators P and N from the design matrix ${\boldsymbol{ \mathcal A }}$ such that

$$\begin{eqnarray}&&{{\boldsymbol{y}}}_{\bullet }={\boldsymbol{P}}\,{\boldsymbol{y}}\end{eqnarray} \tag{ 3 }$$

is the component of the surface map in the preimage and

$$\begin{eqnarray}&&{{\boldsymbol{y}}}_{\circ }={\boldsymbol{N}}\,{\boldsymbol{y}}\end{eqnarray} \tag{ 4 }$$

is the component of the map in the null space. The preimage operator P transforms a vector y in the surface map basis in such a way that it preserves the information in y that gets mapped onto the light curve f via ${\boldsymbol{ \mathcal A }}$ (the preimage) and discards the rest. The null-space operator N does the opposite: it preserves only the information in y that gets mapped onto the zero vector via ${\boldsymbol{ \mathcal A }}$ (the null space). In other words, the P and N operators reveal the components of the surface map that contribute to the light curve ( y _•) and the components that do not ( y _◦). The vector y _• represents all of the information that can be learned from a stellar light curve, while y _◦ represents all of the information that cannot. Importantly, while we carry out this decomposition in the spherical harmonic basis, the null space and preimage are independent of the choice of basis; see the Appendix for a discussion.

It is instructive to visualize these components in an actual surface mapping exercise. Figure 4 shows the decomposition of four hypothetical surfaces (left) into preimage (middle) and null-space (right) components under the flux operator ${\boldsymbol{ \mathcal A }}$, which we compute for definiteness at an inclination of I = 60° over a full rotation. Surfaces are shown in an equal-area Mollweide projection alongside the corresponding light curves. Note that both the true map and its associated light curve are simply equal to the sum of the preimage and null-space components (Appendix). As expected, all of the information in the light curve comes from the preimage, and the null space contributes exactly zero flux at all phases. However, most of the information about the surface is stuck in the null space.

Zoom In Zoom Out Reset image size

In the first row of the figure, corresponding to a surface with a single large spot, it is clear that the light curve contains information about the presence of the spot at roughly the correct longitude and latitude. There are additional artifacts across the stellar surface, and the exact shape of the spot is not well constrained by the data; these issues, however, can easily be resolved with a circular spot prior.

The degeneracies of the mapping problem are much more apparent in the second and third rows, corresponding to surfaces with many smaller spots. The locations, sizes, and very existence of most of the spots are simply not encoded in the light curve. Even with an extremely restrictive prior, it may be difficult—if not impossible—to learn the properties of the spots from the individual light curves.

The fourth row corresponds to a surface with much higher resolution features; this example highlights how the information content of light curves all but vanishes at small scales. Virtually all of the spatial information at the scales of interest is in the null space.

The model mapping spherical harmonic coefficients to observed fluxes described in the preceding sections is linear, so we might expect that we could fit for the spherical harmonics using linear least-squares,

$$\begin{eqnarray}&&\hat{{\boldsymbol{y}}}\mathop{=}\limits^{?}{\left({{\boldsymbol{ \mathcal A }}}^{\top }\,{\boldsymbol{ \mathcal A }}\right)}^{-1}\,{{\boldsymbol{ \mathcal A }}}^{\top }\left({\boldsymbol{f}}-{\bf{1}}\right).\end{eqnarray} \tag{ 5 }$$

But, since the flux operator ${\boldsymbol{ \mathcal A }}$ is low-rank, the model is underdetermined, so the above operation is not defined, and, equivalently, the Fisher information matrix

$$\begin{eqnarray}&&{\boldsymbol{ \mathcal I }}({\boldsymbol{y}})=\displaystyle \frac{1}{{\sigma }_{f}^{2}}\,{{\boldsymbol{ \mathcal A }}}^{\top }\,{\boldsymbol{ \mathcal A }}\end{eqnarray} \tag{ 6 }$$

is singular. In Equation (6), σ_f is the flux measurement uncertainty. There exist many methods and procedures for linear regression for underdetermined models, but all solutions amount to imposing stronger assumptions about the model in order to break the degeneracies.

One standard method is regularized least-squares, in which Equation (5) becomes

$$\begin{eqnarray}&&{\hat{{\boldsymbol{y}}}}_{\lambda }\equiv {\left({{\boldsymbol{ \mathcal A }}}^{\top }\,{\boldsymbol{ \mathcal A }}+\lambda \,{\boldsymbol{I}}\right)}^{-1}\,{{\boldsymbol{ \mathcal A }}}^{\top }\left({\boldsymbol{f}}-{\bf{1}}\right),\end{eqnarray} \tag{ 7 }$$

where I is the identity matrix and λ is a parameter controlling the strength of the regularization. In a Bayesian context, this can be interpreted as placing a Gaussian prior with variance ${\sigma }_{0}^{2}=\displaystyle \frac{{\sigma }_{f}^{2}}{\lambda }$ on the spherical harmonic coefficients.

In the limit of infinitesimal regularization, it can be demonstrated (see, e.g., Hogg & Villar 2021) that ✓

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\lambda \to {0}_{+}}{\hat{{\boldsymbol{y}}}}_{\lambda }={{\boldsymbol{y}}}_{\bullet }.\,\end{eqnarray} \tag{ 8 }$$

In this sense, y _• represents our knowledge about the surface of a star after an observation if we have no prior information whatsoever on y . Any other information about the surface is entirely driven by our assumptions.

Equation (8) is the solution to the surface map we would obtain if we knew nothing about the stellar surface before analyzing the light curve. In practice, this is never really the case. For instance, we know that stellar surfaces must have nonnegative intensities everywhere. While this may seem like a trivial prior, nonnegativity can be a powerful degeneracy-breaking constraint (e.g., Fienup 1982). Moreover, we know that stellar surfaces usually consist of localized features; under an appropriate compactness prior, solutions like the preimage in the fourth row of Figure 4 (for example) could be confidently ruled out.

The effect of different modeling assumptions on the structure and size of the null space is beyond the scope of this paper. In general, this is a particularly difficult question to address because non-Gaussian priors on the surface map break the linearity of the problem. While compactness constraints (like the assumption of a small number of discrete circular spots) can break many of the degeneracies discussed here, simulation-based arguments show that it is still not possible to uniquely constrain their number, locations, or sizes from individual light curves (Basri & Shah 2020).

Figure 5 shows the same decomposition of stellar surfaces into what can (the preimage) and cannot (the null space) be learned from a light curve, but this time for stars viewed at an inclination I = 85°. Interestingly, the structure of the null space is somewhat different; in particular, features in the preimage are latitudinally smeared. This issue is well known (e.g., Cowan et al. 2009; Basri & Shah 2020), and any information about which hemisphere a feature is in formally vanishes as I → 90°.

Zoom In Zoom Out Reset image size

Instead of lamenting the difficulties of constraining stellar surface features at near-edge-on orientations, let us focus on the fact that the null space is a function of inclination. A useful property of the spherical harmonics is that a rotation operation on any component in the basis can only change the order m of the harmonic; the degree l is constant under rotation. In other words, rotating any of the spherical harmonics in Figure 2 about an arbitrary axis simply yields a weighted sum of the spherical harmonics along its row. ¹⁴ This means that changing the inclination of the star—which changes the axis of rotation in the observer's frame—simply changes the weighting of modes that give rise to certain signals in the light curve. This, in turn, results in the dependence of the null space on inclination.

To better understand this effect, let us define the quantity

$$\begin{eqnarray}&&S\equiv 1-\mathop{\mathrm{lim}}\limits_{\lambda \to {0}_{+}}\mathrm{diag}\left(\lambda {\left({{\boldsymbol{ \mathcal A }}}^{\top }\,{\boldsymbol{ \mathcal A }}+\lambda \,{\boldsymbol{I}}\right)}^{-1}\right),\end{eqnarray} \tag{ 9 }$$

which we will refer to as the variance reduction of the coefficients y characterizing the stellar surface. Note that in a Bayesian context, this is equivalent to the posterior shrinkage

$$\begin{eqnarray}&&S=1-\mathop{\mathrm{lim}}\limits_{{\sigma }_{0}^{2}\to \infty }\displaystyle \frac{{\sigma }^{2}}{{\sigma }_{0}^{2}},\end{eqnarray} \tag{ 10 }$$

where ${\sigma }_{0}^{2}$ is the prior variance and σ² is the posterior variance on a particular surface mode we are trying to constrain (see, e.g., Betancourt 2018).

The variance reduction is a dimensionless quantity describing how informative a measurement is about a given mode on the surface in the limit of an infinite signal-to-noise ratio (S/N) and is independent of what the stellar surface actually looks like. If, at an infinite S/N and with no regularization (or, in a Bayesian context, a completely uninformative prior), a particular mode can be learned exactly from a data set, the variance reduction is defined to be unity. Conversely, if the data are completely unconstraining of that mode (i.e., it is entirely in the null space), S will tend to zero.

Figure 6 shows the variance reduction S given a single observation of a star at a random inclination. Each thin blue curve corresponds to a particular draw from an isotropic inclination distribution; the thick blue curve is the average over 300 trials. For some of the low-degree modes, S is relatively high; it is fairly easy to constrain the dipole moment from a light curve, as this is usually the dominant sinusoidal signal. However, as the degree l increases, S decreases dramatically; at l = 14, corresponding to features on scales of roughly 13°, the light curve can only tell us ∼10% of the total information about what the surface looks like. As l increases further, S tends to zero. Another important feature of S, which we hinted at above, is that it is exactly zero for all odd-degree modes above l = 1. This is a well-known fact; all odd spherical harmonics other than the dipole are in the null space regardless of inclination (e.g., Luger et al. 2019). In other words, these spherical harmonics are perfectly antisymmetric in projection over the unit disk when viewed from any orientation. Absent structure to break these symmetries (see below), we simply cannot learn anything about these modes from stellar light curves. If we average over S for all modes up to l = 15, we find that a single light-curve measurement can only tell us ∼9% of the information about the surface on those scales.

Zoom In Zoom Out Reset image size

Fortunately, however, there is quite a bit of scatter in S for different values of the stellar inclination. As we will see, we can use this dependence of the null space on inclination to our advantage. If we could observe a star from many different vantage points, we would be able to break many of the degeneracies at play, since we would get different constraints on the amplitude of each mode when viewed at different inclinations. This, of course, is not possible (at least not yet). But what we can do is observe many similar stars, each viewed at a different (random) inclination, and attempt to learn something about the properties of the ensemble of stars as a whole. In the following section, we explore the role of ensemble analyses in breaking the degeneracies of the mapping problem in more detail.

In an ensemble analysis, we assume that we observe the light curves of many stars that are "similar" in some statistical sense. As a thought experiment, let us consider an extreme version of ensemble analysis in which all of the stars in our sample happen to have identical surfaces. We will still assume they are oriented at random inclinations, as we would expect for field stars. Figure 7 shows the variance reduction curves for this hypothetical scenario, assuming we have access to light curves of 1 (blue), 3 (orange), 10 (green), and 30 (red) identical stars viewed at random inclinations.

Zoom In Zoom Out Reset image size

The addition of light-curve measurements at different orientations increases the variance reduction at all even spherical harmonic degrees (the odd degrees, as we mentioned above, are always invisible). Note that, since we are in the limit of infinite S/N, the fact that we have more light curves (i.e., more data) is irrelevant; the increase in the variance reduction is instead due to the fact that our observations from different vantage points broke some degeneracies in the problem. This is a consequence of the fact we mentioned in the previous section: the null space (for the even modes) is a strong function of the inclination.

If we average over all modes, we obtain an average S of ∼24% for l ≤ 15 when our sample size is three (orange curves); we have more than doubled the information content of our observations. If we further increase our sample size to 10 (green curves), the variance reduction approaches 100% for all even l ≤ 15 modes, effectively saturating for a sample size of 30 (red curves). Thus, if we were able to measure the light curves of identical stars from many different inclinations, the null space would consist only of the odd modes. In the limit of a large number of light curves, and assuming all stars in the sample have identical surfaces, ensemble analyses can tell us up to half of all of the information about the stellar surfaces.

To understand what we can learn about the surfaces in this limit, let us return to the stellar surfaces we considered in Figures 4 and 5. Figure 8 shows the same decomposition of these surfaces into preimage and null space, but this time assuming we measure the light curves of these stars from many different random inclinations. As expected, the partition of information between the preimage and the null space is about 50–50. Because spots are compact features, they are necessarily made up of a continuum of spherical harmonic modes spanning many different values of l; they can therefore be seen in both the even (the preimage) and odd (the null space) modes. The absence of information about the odd modes therefore does not affect our ability to infer the shape and location of the features on the surface. ¹⁵ Interestingly, however, the symmetries at play require spots to be paired with antipodal dark mirror images in the preimage and with bright ones in the null space (which sum to perfectly cancel out in the true map). Thus, there are still degeneracies in this very idealized ensemble problem, but they are much easier to break with a suitable choice of prior. For instance, in the single-spot case (first row), the "ghost" image in the southern hemisphere is surrounded by a bright ring (whose effect is to cancel out its contribution to the flux); either a compactness prior or a prior that enforces uniformity in the background could easily penalize that feature in the fit. This may be much harder to do for surfaces like that shown in the second row, but we can still (in principle) learn about the size, shape, and latitude (if not the number) of starspots from the light curves.

Zoom In Zoom Out Reset image size

There are several practical reasons why this kind of ensemble analysis may be very difficult (but not impossible) in practice. Before we address those, let us consider an important effect that we have not addressed thus far in this paper: limb darkening.

Typically, the shallower the angle between the line of sight and the stellar surface normal, the higher up in the stellar atmosphere optical depth unity is reached. At optical wavelengths, lines of sight directed toward the limb of the star therefore probe cooler temperatures, resulting in the well-known effect of limb darkening. Features close to the limb of the star therefore contribute less to the total outgoing flux, and this must be accounted for when computing the effect of a rotating starspot on the light curve. Limb darkening is often parameterized as a low-order polynomial in the cosine of the line-of-sight angle (Kopal 1950).

To understand how limb darkening affects our ability to infer surface properties from stellar light curves, let us repeat the experiment from the previous section, this time with moderate quadratic limb darkening (with coefficients u₁ = 0.5 and u₂ = 0.25, although our conclusions do not qualitatively change for different values). Figure 9 shows the variance reduction plot (same as Figure 7 but this time accounting for limb darkening). Interestingly, there is no longer a clean division of the null space between even and odd modes in the limit of a large number of light curves. This is because limb darkening effectively lifts the odd modes out of the null space at the expense of the even modes. While no coefficient lies entirely in the null space (S = 0) when limb darkening is present, no coefficient can be uniquely inferred (S = 1) either. This can be understood by noting that a polynomial limb-darkening law can be written exactly as a linear combination of the m = 0 spherical harmonics up to a degree equal to the order of the limb darkening (Luger et al. 2019; Agol et al. 2020; in this case, l = 2). Since the limb-darkening operation is a (multiplicative) downweighting of the surface intensity, the map seen by the observer is just the product of the spherical harmonic representations of the surface ( y ) and the limb-darkening profile. And since spherical harmonics are just polynomials on the surface of the sphere, the product of the spherical harmonics of degree l₁ and l₂ is a spherical harmonic of degree l₁ + l₂. This means that the linear limb-darkening component (l = 1) effectively raises the degree of all spherical harmonic coefficients of the surface map by 1. This has the effect of reversing the null space; under only linear limb darkening, it is the even modes that would be in the null space. However, the quadratic limb-darkening term (l = 2) raises the degree of all spherical harmonics by 2, so its presence ensures that the even modes can still be probed to some extent. In reality, the true limb-darkening profile of a stellar surface is more complicated than a two-parameter quadratic model can capture; but one may still expand it as an arbitrary order polynomial, in which case the argument still applies—limb darkening mixes the null space and the preimage in a nontrivial way. The fact that no coefficient can be determined uniquely—i.e., there are perfect degeneracies involving all modes on the surface—could make it more difficult in practice to perform ensemble analyses on limb-darkened stars.

Zoom In Zoom Out Reset image size

In reality, it is unlikely that all stars in a given ensemble will have exactly the same limb-darkening coefficients, however "similar" the stars may be. Figure 10 shows the same variance reduction plot as Figure 9 but for the case where each star has coefficients u₁ = 0.5 ± 0.05 and u₂ = 0.25 ± 0.025; i.e., we add a scatter of 10% in the value of these coefficients. The plot shows the variance reduction in the hypothetical case where we know the coefficients for each star exactly. Now, as the size of the ensemble increases, the variance reduction approaches unity for all spherical harmonic modes. In the same way that the null space is a strong function of the inclination, allowing us to chip away at it with observations at different inclinations, the null space is also a strong function of the limb-darkening law. Even a small amount of variance in the coefficients is sufficient to constrain all surface modes exactly (in the limit of infinite S/N and a large number of light curves).

Zoom In Zoom Out Reset image size

In practice, of course, we will never know the limb-darkening coefficients exactly. In the next section, we will revisit this and other assumptions we made above in a more sober light.

There are three major points that make the kind of ensemble analyses discussed above difficult in practice. First, and perhaps most obviously, stars are not identical, no matter how "similar" we think they may be. Even stars of the same spectral type, age, and metallicity will in general have different configurations of spots on their surfaces. When we perform an ensemble analysis on the light curves of a heterogeneous group of stars, we are learning something about the distribution of surface properties across all of the stars in the sample, not the surface properties of any individual star. What exactly we can learn in this case is not immediately obvious and requires a detailed investigation. This is the subject of the next paper in this series, Paper II, where we show that, armed with a good model, we can learn a lot about the distribution of the starspot properties of a heterogenous group of stellar light curves.

The second point is that while observing many similar stars at different inclinations can greatly help us learn about their surfaces from a statistical standpoint, our analysis above assumed we knew what the values of the individual inclinations were. In practice, this will usually not be the case. While we may have good priors for some stars (from spectroscopic $v\sin I$ measurements or the assumption that transiting exoplanet hosts are likely to have inclinations close to I = 90°), for the vast majority of field stars, we will not know much a priori. Since the inclination is typically degenerate with the spot latitude (e.g., Walkowicz et al. 2013), this decreases the constraining power of ensemble analyses. However, as we show in Paper II, there is still enough information in ensembles of ≳50 light curves to independently constrain the spot latitudes and the individual stellar inclinations.

The final point concerns limb darkening, which also has a strong effect on the structure of the null space. While limb darkening can help us in the same way as the inclination, in practice, it is likely to be more of a problem, since the use of incorrect limb-darkening coefficients can lead to bias in the spot properties when doing inference. It is therefore extremely important to use reliable limb-darkening models when doing ensemble analyses; we also explore this in Paper II.

To understand the normalization degeneracy, consider how we might go about simulating a stellar surface. We might add a dark spot somewhere on the surface, either by expanding it in spherical harmonics or by gridding up the stellar surface and setting the intensity of pixels within the spot to a low value. To compute the light curve, we integrate over the projected disk at each point in time. The resulting light curve will have strange units, so we might then divide by the integral of the background intensity over the unit disk, so that we are now in what we will call fractional units: the flux as a fraction of the flux we would measure if the star had no spots on it.

The top panel of Figure 11 shows two mock light curves we might compute following the procedure above. The solid blue curve corresponds to the light curve of a star with a single large equatorial spot of contrast c viewed at an inclination I = 60°. The dashed orange curve corresponds to a star with a spot at the same location but half the contrast, plus a large polar spot of comparable contrast. Because the equatorial spot on this star has half the contrast of that on the first star, the peak-to-trough amplitude of the orange light curve is half that of the blue light curve. Moreover, since the polar spot is always in view on this star, the peak flux is itself only about half that of the first star. If we were given these two light curves in these fractional units, we might be able to infer these basic differences between the two stars (setting aside for the moment all the issues with the null space discussed in the previous section).

Zoom In Zoom Out Reset image size

However, we do not observe stellar light curves in fractional units. Instead, we typically observe in units of counts on the detector, which depend (among other things) on various properties of the telescope. But even if we knew all of these things, as well as the true luminosity of the star and the precise distance to it, and we could truly perform absolute photometry, that still would not be enough to correctly calibrate the light curve into the units we need. To convert to fractional units, we would need to know the brightness we would measure if the star had no spots on it. This depends on the brightness of the unspotted photosphere, which cannot be measured directly unless we actually resolve the star (i.e., interferometrically).

Usually, this is not much of an issue. Astronomers typically circumvent this by self-normalizing the data, i.e., dividing the flux by the mean, median, maximum, or some similar statistic of the light curve. This operation folds the unknowability of the true units under the rug and transforms the light curve into a relative measurement of the star's temporal variability. While relative measurements are typically what we are interested in anyway, this normalization procedure can sometimes fool us into thinking we have access to information that is simply not observable in single-band photometry. To understand why, consider the bottom panel of Figure 11, which shows the same two light curves in what we will call relative units. To get the light curves in these units, we followed the common procedure of dividing each by the observed "continuum" level (the maximum flux in the light curve), subtracting unity, and multiplying by 1000, yielding relative fluxes in units of parts per thousand.

The two light curves, which were distinct in the fractional units we used to generate them, are indistinguishable in the relative units in which we observe them. There is absolutely no information in the relative light curves that can differentiate between the two stellar surfaces shown in the figure. The depth of each of the dips cannot tell us about either the contrast of individual spots or the total number of spots; in fact, in single-band photometry, these two quantities are perfectly degenerate with each other.

To explore this point in a bit more detail, let us consider the normalization degeneracy in the context of ensemble analyses. Even though we cannot uniquely infer contrasts or numbers of spots from individual light curves, perhaps we could harness the power of the ensemble. Let us therefore go back to our thought experiment in which we added spots to a stellar surface. Assuming, for simplicity, that all spots have the same contrast c, every time we add a spot, the flux (in fractional units) decreases by an amount proportional to c; so, to first order, we can approximate an individual light curve as

$$\begin{eqnarray}&&{\boldsymbol{f}}(c,n)={\bf{1}}-c\sum _{i=0}^{n-1}{\boldsymbol{g}}({{\boldsymbol{\theta }}}_{i}),\end{eqnarray} \tag{ 11 }$$

where g ( θ _i ) is some (complicated) function of the properties of the ith spot, as well as the stellar inclination, rotation period, etc., which we denote by θ _i . Now consider many stellar light curves drawn from some distribution controlling the stellar and starspot properties with probability density p( θ ). The mean of the distribution of the light curves (still in fractional units) is then given by

$$\begin{eqnarray}\begin{array}{rcl}\mu (c,n) & = & {\rm{E}}\left[{\boldsymbol{f}}(c,n)\right]\\ & = & 1-c\,{\rm{E}}\left[\sum _{i=0}^{n-1}{\boldsymbol{g}}({{\boldsymbol{\theta }}}_{i})\right]\\ & = & 1-c\,n\,\alpha ,\end{array}\end{eqnarray} \tag{ 12 }$$

where ${\rm{E}}\left[\cdots \right]$ denotes the expected value and

$$\begin{eqnarray}&&\alpha \equiv \int {\boldsymbol{g}}({\boldsymbol{\theta }})p({\boldsymbol{\theta }}){\rm{d}}{\boldsymbol{\theta }}\end{eqnarray} \tag{ 13 }$$

is the expected value of g . Similarly, the variance of the distribution may be computed as

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }^{2}(c,n) & = & \mathrm{Var}\left[{\boldsymbol{f}}(c,n)\right]\\ & = & \mathrm{Var}\left[c\,\sum _{i=0}^{n-1}{\boldsymbol{g}}({{\boldsymbol{\theta }}}_{i})\right]\\ & = & {c}^{2}\,\mathrm{Var}\left[\sum _{i=0}^{n-1}{\boldsymbol{g}}({{\boldsymbol{\theta }}}_{i})\right]\\ & = & {c}^{2}\sum _{i=0}^{n-1}\mathrm{Var}\left[{\boldsymbol{g}}({{\boldsymbol{\theta }}}_{i})\right]\\ & = & {c}^{2}\,n\,({\beta }^{2}-{\alpha }^{2}),\end{array}\end{eqnarray} \tag{ 14 }$$

where $\mathrm{Var}\left[\cdots \right]$ denotes the variance,

$$\begin{eqnarray}&&{\beta }^{2}\equiv \int {{\boldsymbol{g}}}^{2}({\boldsymbol{\theta }})p({\boldsymbol{\theta }}){\rm{d}}{\boldsymbol{\theta }}\end{eqnarray} \tag{ 15 }$$

is the expected value of g ², and we used the fact that the variance of the sum of independent random variables is equal to the sum of their variances.

To summarize, the mean and variance of the ensemble of stellar light curves in fractional units is

$$\begin{eqnarray}\begin{array}{rcl}\mu & = & 1-c\,n\,\alpha \\ {\sigma }^{2} & = & {c}^{2}\,n\,({\beta }^{2}-{\alpha }^{2})\end{array}\end{eqnarray} \tag{ 16 }$$

for some complicated functions α and β of the distribution of stellar inclinations, rotation periods, and starspot properties. If our observations were collected in these fractional units, we could uniquely infer the spot contrast c and the number of spots n, since these enter as c n and c² n in the expressions for the mean and variance, respectively, and these are straightforward statistics to compute from the ensemble. However, because the observations are made in relative units, in which we typically normalize to the mean, the amplitudes of the features in the light curves can only tell us about the ratio

$$\begin{eqnarray}&&\displaystyle \frac{\sigma }{\mu }\propto \displaystyle \frac{c\sqrt{n}}{1-{cn}\alpha }\end{eqnarray} \tag{ 17 }$$

for some value of α. In other words, even photometric ensemble analyses may not be able to tell us about the values of the contrast and the number of spots independently.

A direct consequence of this normalization degeneracy is that it may not be possible to uniquely constrain the total spot coverage of a star from single-band photometry without strong prior assumptions. The total spot coverage f_S is simply the (average) area of a spot times the total number of spots divided by the total area of the sphere, which may be expressed as ✓

$$\begin{eqnarray}&&{f}_{S}=\displaystyle \frac{1}{2}\left(1-\left\langle \cos r\right\rangle \right)n,\,\end{eqnarray} \tag{ 18 }$$

given an angular spot radius r. While r may be uniquely constrained from the covariance structure of the data (Paper II), n cannot.

The arguments above are heuristic and based on only the first two moments of the distribution of light curves in an ensemble. It is possible, at least in principle, that higher-order moments of the data could encode information to break the c − n degeneracy, but these are generally more difficult to constrain from the data. It is also possible that we could place a lower limit on the spot contrast based on the normalized light curve. In Equation (11), we implicitly assumed that spots are allowed to overlap; under this assumption, it is possible to double the contrast of a spot by simply adding another spot on top of it. We could therefore generate light curves with arbitrarily large dips by choosing a small value for c and a very large value for n. In reality, spots do not behave in this way; many overlapping spots likely still have the same effective temperature. Therefore, assuming there are no bright spots, the maximum depth of a feature in the light curve (even in relative units) could be used to place a lower limit on the spot contrast. Moreover, in Paper II, we empirically show that, depending on the choice of prior, ensemble analyses may be able to place stronger constraints on the contrast c than on the number of spots n.

Nevertheless, even with these caveats in mind, the fact remains that degeneracies like the polar spot effect are fundamental; recall Figure 11, in which the two stars have identical relative light curves but very different spot coverage fractions. In general, any azimuthally symmetric mode on the surface is in the null space of the normalized light-curve problem. These components contribute a constant value to the flux at all phases, which effectively gets normalized away when we measure the light curve. A polar spot is just one manifestation of this degeneracy. Bands, or band-like arrangements of spots, will be at least partly in the null space of normalized light curves. These features change the total "spot coverage" of the surface but do not affect the light curve in a uniquely measurable way.

There is one final subtle point concerning the normalization degeneracy that merits discussion. The common procedure of normalizing light curves to their mean, median, or maximum level does not only change the units of the data, it changes the very covariance structure of the light curves.

To understand why, let us consider the procedure of normalizing a light curve to its mean value. Whenever we scale our data, we must always be sure to scale the error bars accordingly. Since in this case, we are dividing the flux by the mean, one might imagine that we could simply divide each of the measurement uncertainties by the same amount. However, this is technically incorrect.

Consider the example in Figure 12. The top panel shows 1000 samples from a sinusoid with random phases and a period equal to 10 times the observational window. In this limit, each light curve is approximately linear, which causes its mean value to roughly coincide with the midpoint of the observation window. Division by the mean value (bottom panel) results in points near the midpoint being driven to unity and points near the edges (whose values differ the most from the mean) to be driven to both large and small values. If our error bars in the original data were uniform (homoscedastic), the error bars in the normalized light curves are not; the standard deviation (or variance) of the data is now distinctly dependent on the phase.

Zoom In Zoom Out Reset image size

While the example shown in the figure is fairly extreme, the idea here is quite general: the normalization procedure changes the covariance structure of the data. In most cases, the nonstationarity (i.e., the phase dependence) of the variance will be quite small. The effect is primarily important for light curves with periods much longer than the observation window. In these cases, not accounting for this effect could introduce bias in light-curve analyses. A detailed investigation of this effect is beyond the scope of this paper, but we do present a method to correct the covariance matrix of normalized light curves for this issue in Paper II.

In this paper, we explored various degeneracies in the problem of inferring a stellar surface map from its rotational light curve. We discussed the idea behind the null space, the set of surface modes that have exactly zero effect on the observed light curve (Section 2). For rotational light curves, we showed that the vast majority of the information about the surface intensity is in the null space and therefore cannot be inferred from unresolved photometric measurements. We showed, in particular, that the size of the null space grows quadratically with increasing spatial resolution, while the number of independent degrees of freedom in the light curve only grows linearly. Consequently, the information content of the light curves is small for large-scale surface features and vanishingly small for small-scale features.

A direct consequence of this fact is that results based on analyses of individual light curves are often extremely sensitive to the particular assumptions of one's model. We therefore urge lots of care in these kinds of analyses, where assumptions of uniform-contrast, discrete, and circular spots or, alternatively, of maximum entropy or "simplicity," may result in significant bias. In this respect, we agree with the recent paper by Basri & Shah (2020), which cautions against the association of individual dips in light curves with individual, discrete spots.

One of the main results in this paper concerning the null space is its dependence on stellar inclination (Section 2.4). Because the modes that lie in the null space change depending on the viewing angle, observations of a star at many inclinations would break many of the degeneracies in the mapping problem. While this is obviously infeasible in practice, the dependence of the null space on inclination motivates ensemble analyses of many similar stars as a way of circumventing the mapping degeneracies by providing a strong data-driven prior. As we show in Paper II in this series, the joint analysis of tens to hundreds of light curves can uniquely constrain the distribution of starspot sizes and latitudes among the stars in the sample. In that paper, we show that ensemble analyses can even break the latitude–inclination degeneracy (e.g., Walkowicz et al. 2013), allowing one to infer individual stellar inclinations, typically to within about 10°.

We also showed how the null space is a strong function of the stellar limb darkening. While this can again be used to our advantage—by harnessing the variance in the limb-darkening parameters across an ensemble of stellar light curves—in practice, it is likely to complicate inference, since any bias in the assumed limb-darkening coefficients will likely result in bias in the inferred surface properties. We revisit this point in Paper II.

One final point that we did not address thus far concerns occultations. Our results regarding the null space apply strictly to rotational light curves, in which all points on the projected disk contribute to the measured flux. During an occultation by a transiting planet or a binary companion, the weighting of surface modes giving rise to the light curve changes substantially. In fact, the presence of an occultor breaks many of the perfect symmetries that give rise to a null space in the first place (e.g., Luger et al. 2019). This fact can be used to infer the properties of spots in the path of the occultor, as was done (for example) by Morris et al. (2017). A detailed investigation of the null space for the occultation problem is beyond the scope of this work.

The second major source of degeneracies is the fundamental unknowability of the true normalization in single-band photometry (Section 3), which is summarized in Figure 11. In a nutshell, the relationship between the amplitude of a feature in the light curve and the contrast of the feature that gave rise to it depends on quantities like the unspotted photospheric brightness, which is not an observable in single-band photometry. This leads to the possibility of distinctly different stellar surfaces having identical relative light curves, as demonstrated in the figure.

In practice, this degeneracy usually manifests as a nonlinear correlation between the contrast c of a spot and the total number n of spots on the surface of the star. Even in ensemble analyses, one cannot generally learn these two quantities independently from the data alone, only a (complicated) function of the two (Equation (17)). In Paper II, we show empirically that careful analysis of the covariance structure of ensembles of light curves may shed some light on c but cannot uniquely constrain n. A direct consequence of this degeneracy is that constraints on the total number of spots or the spot coverage of a star from single-band light curves depend strongly on the model assumptions.

Recently, Morris (2020b) performed an ensemble analysis of Kepler, K2, and TESS light curves to derive a relationship between the fractional spot coverage f_S and stellar age. That work found that the spot coverage f_S as a function of age is well modeled by a simple power law, decreasing from ∼10% for the youngest (∼10 Myr) stars to less than 1% for the oldest (∼5 Gyr) stars. While this broadly agrees with the expectation that stellar magnetic activity decreases over time, our work strongly suggests that these results depend heavily on the prior. This is because the expression for f_S (Equation (18)) depends on two quantities: the average spot radius r, which can be constrained (see Luger et al. 2021), and the total number of spots n, which we showed cannot be uniquely constrained from single-band light curves. In fact, Morris (2020b) assumed n = 3 for simplicity when doing posterior inference. We therefore urge care in interpreting those results, as it is possible that this assumption does not hold across the large range of spectral types and stellar ages considered in that study. Nevertheless, the other central result of that paper, the relationship between the "smoothed amplitude" of a light curve and the stellar age, is valid. Moreover, the core idea in Morris (2020b; and related studies such as Jackson & Jeffries 2013) is very similar to that advocated here: the use of ensemble analyses to constrain population-level parameters when individual data sets are not sufficiently constraining.

Another recent paper relevant to our work is that of Basri & Shah (2020), who investigated the information content of stellar light curves, exploring what can and cannot be learned about starspot configurations from individual light curves. That paper strongly urges against the common practice of interpreting light curves with one or two dips as originating from one or two spots, respectively, a point we strongly agree with (see our Figure 1). It also reinforces our point about the additional degeneracies introduced by the unknown normalization inherent to single-band photometry.

As with the degeneracies due to the null space, there are potential ways to break the normalization degeneracy. In principle, the maximum level of a light curve could set a lower limit on the brightness of the unspotted photosphere, particularly for long-baseline, time-variable surfaces (Basri 2018). However, this will work only if the surface is known to be made up exclusively of dark spots. The presence of bright spots (faculae), which are common on the Sun, make it difficult for one to infer this quantity (and hence the correct normalization) from single-band photometry in practice.

A better approach may be to collect photometric data in multiple wavelength bands, an idea that has been explored recently (e.g., Gully-Santiago et al. 2017; Guo et al. 2018). Assuming the locations and sizes of a star's spots are constant in wavelength, the amplitude of the light curve in different bands (and, in particular, its slope as a function of wavelength) can be used to directly constrain the temperature, and hence the contrast, of the spots. This effectively breaks the c − n degeneracy. In practice, the effective sizes of spots may be different at different wavelengths, which could complicate this picture somewhat, but the point still stands that light curves collected in multiple bands contain at least partial information about the correct normalization. A more detailed analysis of the information content of multiband light curves and, by extension, spectroscopic time series is deferred to a future paper in this series.

There are several simplifying assumptions we made in this paper that are worth discussing. First, in our characterization of the null space, we assumed we knew the stellar inclination, stellar rotation period, and limb-darkening coefficients exactly, and we computed the information content of stellar light curves in the limit of an infinite S/N. None of these assumptions are valid in practice. In realistic scenarios, the variance reduction will be necessarily lower than what we presented here; the curves in Figure 7, for instance, are therefore strict upper bounds on the amount of information that can be learned about the surface. In some cases, it may be possible to infer the inclination from spectroscopic $v\sin I$ measurements or the limb-darkening coefficients from the shape of transiting exoplanet light curves, but any uncertainty in these quantities will degrade our ability to constrain information about the surface.

Moreover, our constraints on the information content of ensembles of light curves are also strict upper limits, since we assumed that all stars in the ensemble are identical. Variance in the surface maps of stars in a population makes it effectively impossible to infer detailed properties of the individual stellar surfaces in the ensemble. We argued without proof that what we can learn instead are the properties of the distribution of stellar surfaces among the population, such as the distribution of spot sizes and active latitudes. Demonstrating this point is more difficult, so we defer it to Paper II, where we construct a custom GP model and show that it can be used to infer such population-level parameters from ensembles of light curves.

Finally, we limited our discussion to static stellar surfaces. Real stellar surfaces vary with time as individual spots form, evolve, and dissipate or a star progresses through its activity cycle. This, too, makes it more difficult to constrain the stellar surface at a particular point in time, since the amount of data corresponding to a particular surface configuration is more limited. On the other hand, a time-variable surface could be used to our advantage from an ensemble standpoint. Although we are unable to uniquely infer individual spot properties, we could treat each point in the light curve as a realization of a random process (i.e., a draw from some distribution) describing the stellar surface at a high level. In this sense, we could harness the fact that we can measure the light curves of many "similar" surfaces to learn something about their statistical properties. Demonstrating this is beyond the scope of the present work, but we return to this point in Paper II.

This paper is the first in a series dedicated to developing methodology to infer stellar surface properties from unresolved measurements. It sets the stage for Paper II, in which we develop an interpretable GP model for starspot-induced light-curve variability. This GP is aimed specifically at the difficult problem of jointly modeling many light curves in an ensemble analysis. Although we hinted at this here, we will show in that paper that ensemble analyses can uniquely constrain several statistical properties of starspots, including their distribution of radii and latitudes across stars in the ensemble. In that work, we also consider temporally evolving surfaces, which we have neglected in our discussion thus far.

Finally, while this paper explicitly deals with the problem of photometric rotational light curves, the degeneracies outlined here and the methodology developed in this series of papers apply more broadly in other contexts. These include applications where the stellar surface is a nuisance (exoplanet detection and characterization using transit light curves or radial velocities) and spectral time series data sets (such as transmission spectroscopy and Doppler imaging).

In keeping with other papers in the starry series, all figures in this paper are generated automatically from open-source scripts linked to in each of the captions (⏬), and the principal equations link to associated unit tests that ensure the accuracy and reproducibility of the algorithm presented here (✓/×).