\section{Introduction}

The inferred high spherical albedo values of several hot super-Earths \citep{Demory2014,Malavolta2018} are well  above what the current literature predicts for magmatic surfaces \citep{Essack2020,Rouan2011,Edgett1997}. Two main explanations have been put forward: reflective exotic magmas \citep{Rouan2011} and reflective volcanic/mineral atmospheres \citep{Hamano2015,Pluriel2019}. However, it is currently not possible to verify which mechanism applies to specific super-Earths due to observational data being either limited or degenerate. Generally, constraining the properties of super-Earth exoplanets is troublesome due to their small sizes \citep{Valencia2013,Dorn2017}. In spite of this, advances in mass-radius measurements and superior spectroscopic analyses of their atmospheres \citep[e.g.][]{Tsiaras2016,Esteves2017} has lowered their compositional degeneracies. These advancements have permitted researchers to apply geological models to specific super-Earths in order to better understand their structures \citep[e.g.][]{Zeng2013,Zeng2016,Zeng2018}. However, one aspect that remains unexplored is the surface composition and topology of these bodies. Directly analysing the geomorphology of extra-solar planets is not possible due to the observational limitations set by the \textit{Rayleigh Criterion}. For example, if one wanted to resolve a feature that occupies $10\%$ of a planetary surface (e.g. an ocean) on an Earth-sized exoplanet located $5~\rm pc$ away, a spatial resolution of a few microarcseconds would be required. However, this is beyond our present capabilities which is one of the reasons why it is troublesome to decipher which mechanisms lead to large albedos on molten super-Earths such as Kepler-10b, Kepler-21b and K2-141b \citep{Rouan2011,Demory2014,Malavolta2018}.

The purpose of our study is to investigate the spherical albedo values of airless and near-airless, molten super-Earths with a particular focus on its dependence to the composition and ocean waviness. As this study is a first stepping-stone to theoretically estimate the albedo of these ultra-hot exoplanets, we had to set some simplifications (see Sec.~\ref{sec:limitations}) and left some tasks for future works (see Sec.~\ref{sec:future_work}). We show that the ocean's roughness will lead to a higher degeneracy in the interpretation of the spherical albedo values of these systems. The close-orbiting super-Earths relevant to our study will be tidally locked therefore containing a molten day-side and a cold night-side. Due to the high day-side temperatures, a vaporised mineral atmosphere will form that rapidly traverses into the night-side where the temperatures are low enough for the minerals to condense \citep{Schaefer2012,Ito2015,Kite2016}. Whilst the atmosphere is expected to be generally optically thin in most of the infrared spectral range, the rapidly expanding gas would carry enough kinetic energy to trigger magma ocean waves (see Sec.~\ref{sec:hurst} of the appendix and \citet{Kite2016}). In Sec.~\ref{sec:choice_of_planet} we justify our choice of using Kepler-10b as an analogue planet for our spherical albedo model. In Sec.~\ref{sec:generate_surfaces} we model the magma ocean waves using fractal mathematics. In Sec.~\ref{sec:atmospheric_effects} we explain that while mineral atmospheres may affect the spherical albedo values within the spectral range of the \textit{Kepler} telescope, for future missions with larger wavelength ranges, the effects would be negligible. In Sec.~\ref{sec:Results} we show the results of our simulations and we present an analytic approximation for the planar albedo, which we then integrate across all angles to determine the spherical albedo. Our analytic functions model the spherical albedo values of molten, airless magma ocean super-Earths for different roughnesses and surface compositions. Finally, in Sec.~\ref{sec:Discussion} and \ref{sec:conclusion} we discuss and present our results respectively.

Note that many of the terms used within this manuscript to describe the reflective properties of super-Earths are similar but not identical so it is necessary to define them to avoid confusion. The reflectance, $R$, is a directional quantity that depends on the incident and emerging angles, as well as the wavelength of the light and the composition (via the refractive indices) and roughness of the surface. The spherical albedo, $A_{S}$, is the total fraction of incident radiance scattered by a body into all direction. The planar (or plane) albedo, $A_{pla}$, is the spherical albedo of a planar surface of the planet. Finally, the bolometric Bond albedo, $A_{B}$, is the spherical albedo integrated across all wavelengths.

\section{Our Choice of Planet}
\label{sec:choice_of_planet}

Throughout the rest of this manuscript we will use Kepler-10b as an analogue for molten, airless super-Earths within our spherical albedo models. We do this for two reasons: first we want to input realistic parameters into our simulations instead of speculative idealised ones. Secondly, we want to show how our analytic function for the spherical albedo values can be applied with ease to any molten, airless super-Earth of choice. Therefore, before we progress into an explanation of our models we would first like to justify our choice of using Kepler-10b as our test planet. We believe this body is an adequate candidate for our method because most mass and radius measurements are consistent with a rocky composition \citep[e.g.][see Fig.~\ref{fig:kepler10b}]{Batalha2011,Fogtmann2014,Dumusque2014,Esteves2015,Weiss2016,Dai2019}.
\begin{figure}[h]
	\centering
	\includegraphics[scale=0.8]{Kepler_10_b_mr.png}
	\caption{The mass and radius of Kepler-10b according to different studies. The planetary composition model curves are from \citet{Zeng2013} and \citet{Zeng2016}. The variations in the measurements are due to a combination of uncertainties in the stellar parameters and using different observational raw data.}
	\centering
	\label{fig:kepler10b}
\end{figure}
Furthermore, the strong temperature contrast inferred from the observed phase curve between the day-side temperature of $\simeq 2750$~K and night-side temperature $\simeq 50$~K implies that a thick atmosphere is improbable \citep{Rouan2011}. In addition, the large thermal gradient is suggestive that Kepler-10b is tidally locked. Hence, the vaporised mineral gases formed on the day-side \citep{Schaefer2012,Ito2015} will rapidly expand into the night-side and condense, forming waves in the magma ocean beneath as they travel \citep{Kite2016}. Whilst there is still an ambiguity in the measured data, in this paper we focus on the scenario where Kepler-10b is a rocky exoplanet hosting a fully magmatic day-side that is lacking a significant atmosphere. Currently, the best fit value for the bolometric Bond albedo is $0.48 \pm 0.35$ \citep{Rouan2011}, which is the one we will reference throughout this paper.
This study investigates the spherical albedo values for various magma ocean compositions and surface roughnesses. The spherical albedo is useful to consider because it is defined to be a monochromatic property, so it is independent of the irradiation spectrum of the host star. From our spherical albedo values and an estimation of the stellar spectra, the bolometric Bond albedo values could be inferred \citep[see][for an example of a similar concept but for planets with atmospheres]{Marley1999}. It is also important to note that the primary objective of this study is to present our model for constraining the composition of exoplanetary magma oceans; we are therefore not focusing on formational or observational limitations.

\section{Models}
\subsection{Generating our Magma Ocean Fractal Surfaces}
\label{sec:generate_surfaces}

One commonly used, although incomplete, definition of a fractal was given by Benoit Mandelbrot as a shape whose fractal dimension (also called Hausdorff dimension) is greater than its topological dimension. The fractal dimension is a measure of the geometric complexity of a system so that from a geological perspective, this parameter would provide information on the roughness of a surface. For instance, Earth’s oceans have calculated fractal dimensions of $\approx2.3$ \citep{Stiassnie1991,Stiassnie1991b} which means that due to surface inhomogeneities Earth’s true surface area could be several orders of magnitude greater than what is predicted from a perfect spherical model. Similar analyses have also been performed on the fractal-nature of Mars \citep{Deliege2017,Demin2017} and Venus \citep{Demin2018}. Roughness has also been shown to affect the reflective properties of non-fluid surfaces with coarser systems having lower ones than their smoother counterparts \citep[e.g.][]{Bennett1961}. This is because with a rough surface photons can get trapped in crevices and consequently absorbed. In contrast, with smooth surfaces photons are absorbed less as the albedo values are only dependent on the material properties. For example, soil \citep{Matthias2000} and ice \citep{Lhermitte2014} have been found to be almost $\sim 50\%$ less reflective than their standard values when they are rough. This is explained by a difference in the phase function of a rough surface and that of a smooth one with the same composition. In addition, several studies have found strong linear relationships between soil's albedo values and its roughness \citep{Oguntunde2006,Cierniewski2013}. In contrast, when dealing with fluids their bolometric Bond albedo values do not scale proportionally to the roughness of the system \citep{Haltrin2001}. This is because liquid waves are topologically distinct from the types of concavities one would find in solids (e.g. cracks, holes etc.). Waves are not as efficient in trapping photons so other effects such as the angle of incidence of the light ray relative to the immediate surface, and the number of internal reflections are important as they determine how much energy is lost. Since this cannot be modelled analytically, in this paper we perform ray tracing simulations in order to determine the reflective properties of synthetic magma ocean surfaces.

To keep the simulations inexpensive we will use Schlick's approximation, Eq.~\ref{eq:Schlicks}, for the reflectance at a liquid boundary \citep{Schlick1994}. Schlick's approximation is not an exact solution, but rather a method of keeping the ray tracing simulations inexpensive, whilst still maintaining physical accuracy. This method only works when the ray of light is traversing from a medium with a low refractive index to that of a higher one, which is relevant to our model. We are aware that this approximation will inevitably result in a decrease of accuracy but we unfortunately had to make this compromise in order to cover the large parameter space within a reasonable time-frame. In Sec.~\ref{sec:comparison_of_models} of the appendix we compare Schlick's approximation with other models and the exact Fresnel equations in order to quantify their differences. We define the parameter $n_{sy}$ as the synthetic refractive index of the magma, which accounts for the effects of the real ($n$) and imaginary ($k$) indices. The derivation of the synthetic refractive index, $n_{sy}$, is shown in Sec.~\ref{sec:derivation_n_sy} of the appendix. Therefore, $n_{sy}$ is dependent on the geochemical composition of the medium and the wavelength of the light as shown by the best-fit equations listed in Table~\ref{tab:bestfitequations}:
\begin{subequations}
	\begin{equation}
	R\left(\theta \right) \approx R_{0} + \left( 1 - R_{0} \right) \left(1 - \cos{\theta} \right)^{5}
	\label{eq:Schlicks}
	\end{equation}
	\begin{equation}
	R_{0} = \left(\dfrac{n_{sy} - 1}{n_{sy} + 1} \right)^{2}
	\label{eq:R0}
	\end{equation}
	\begin{equation}
	\ n_{sy} = \frac{ \left(\frac{ \left|n - 1 \right|^{2} + \left|k \right|^{2} }{ \left|n + 1 \right|^{2} + \left|k \right|^{2} } \right)^{1/2} +1 }{1 - \left(\frac{ \left|n - 1 \right|^{2} + \left|k \right|^{2} }{ \left|n + 1 \right|^{2} + \left|k \right|^{2} } \right)^{1/2} }
	\label{eq:nsynth}
	\end{equation}
\end{subequations}
Where $ \theta $ is the angle of incidence and $R$ is the reflectance. Having chosen our equations, the next step was to model the fractal ocean surfaces. We used the \textit{fbm 0.3.0} python package with the \citet{Davies1987} fractional Gaussian noise (fGn) method. This technique is theoretically exact in generating discretely sampled fGn numerical values for a given fractal dimension. A thorough explanation of the Davies and Harte method can be found in page 412 of \citet{Wood1994}. We chose to model the magma ocean surfaces as a collection of different points with elevations given by fGn in order to approximate the random motions induced from the vapourised atmospheric winds (see Sec.~\ref{sec:hurst} of the appendix for a mathematical relationship of the wind energies and the surface roughness). However, this approach assumes that each point is independent of its neighbours which would only be true for a zero viscosity fluid (i.e. a superfluid) that also lacks surface tension. Clearly magmas have viscosities so the non-continuous, sharp terrains are unrealistic. Fig.~\ref{fig:fractalsurfacesharp} is an example of a pre-smoothed, unrealistic surface.
\begin{figure}[h]
	\centering
	\includegraphics[scale=0.70]{Fractal_terrain_sharp.png}
	\caption{The unsmoothed fractal surface of an ocean whose wave heights are yet to be adjusted for. This surface was generated for a Hurst exponent of 0.7. The circled area is a magnification of a $\times 100$ for the region $0-1~\rm m$}
	\centering
	\label{fig:fractalsurfacesharp}
\end{figure}
In order to smooth these surfaces we applied a Gaussian filter with a standard deviation of 100 pixels (equivalent to 1~m) given by the \textit{Astropy} package \citep{Astropy2013,Astropy2018}. The process of smoothing the surface is analogous to adding a viscosity which, by definition, will act against the fluid flow. We chose to smooth out perturbations smaller than 1~m as this closely matches the waves we observe on Earth's oceans (see Sec.~\ref{sec:limitations} for a justification of why the waves of Earth's oceans are an appropriate analogue for the roughness of the magma ocean on Kepler-10b and other ultra-hot super-Earths). From a mathematical perspective, the act of smoothing perturbations smaller than $1~\rm m$ means that fractal behaviour seizes at very small scales which is expected for natural entities instead of idealised mathematical constructs. It must be noted that maximising the Gaussian filter will produce a smooth surface so that the planar albedo will be given by Eq.~\ref{eq:R0}. If instead no filter is applied then the planar albedo will converge to 1 as the individual waves will be very thin and tall (i.e. the waves will be steep so the angle $\theta$ will approach $90^{\circ}$; see Fig.~\ref{fig:fractalsurfacesharp}). These two scenarios are unrealistic as they require extreme viscosities which is unlikely given the thermodynamic conditions and magma compositions considered in this study. Concerning with the average wave peak amplitude, this required finding a relationship between the Hurst exponent, which measures the coherence of the surface, and the root mean squared (RMS) roughness. Experimental results from \citet{Durst2011} found that the fractal dimension and the RMS were related by the following equation:
\begin{equation}
RMS \approx h_{0}\left(F - 2 \right)
\label{eq:Durst}
\end{equation}
Where $h_{0}$ is a height constant and $F$ is the fractal dimension. For self-similar processes where local and global properties are interlinked, the fractal dimension is related to the Hurst exponent by $F = 3 - H$. Therefore \textit{H} can be substituted in:
\begin{equation}
RMS \approx h_{0}\left(1 - H \right)
\label{eq:DurstH}
\end{equation}
We decided to work with the Hurst exponent ($H$) instead of the fractal dimension ($F$) in order to make our code easier to adopt. This is because the Hurst exponent is independent of the topological dimension, whilst the fractal dimension is not. Therefore, if one wanted to implement our code into a 1-D, 2-D or 3-D radiative transfer model, the same Hurst exponent could be used resulting in the same spherical albedo value for each scenario. The constant $h_{0}$ is calculated experimentally, which can be done using Earth's oceans as an analogue. We know that these have a Hurst exponent of $\approx 0.7$ \citep{Stiassnie1991,Stiassnie1991b} with a root mean squared average waveheight of $\approx 1.5$~m \citep{MetOffice2010}. Consequently, solving Eq.~\ref{eq:DurstH} for $H = 0.7$ and $\rm RMS = 1.5$~m gives $h_{0} \approx 5$~m, yielding the final form of our equation:
\begin{equation}
	RMS \approx 5\left(1 - H \right)
\label{eq:DurstH2}
\end{equation}
In order to correct for the height of our magma ocean waves we implemented Eq.~\ref{eq:DurstH2}. This was done by taking the smoothed (i.e. the Gaussian filter had been applied), unadjusted surface and finding the average (absolute) height. The RMS value (calculated with Eq.~\ref{eq:DurstH2}) was then divided by the average height of the smoothed surface to find the average correction factor. Finally, we multiplied the average correction factor to each point in the smoothed, unadjusted surface in order to get our more realistic magma ocean surface as shown by Fig.~\ref{fig:fractalbouncing}. Note that we are assuming that the parameter $h_{0}$ is constant and equal for the Earth and Kepler-10b. This assumption has been made because the viscosity of magma at very high temperatures is fortuitously similar to that of water. We explain this point more thoroughly in Sec.~\ref{sec:limitations}.

Once we had the fractal surface generator we created a basic simulator that would throw 10,000 light rays whose initial locations were uniformly sampled. Each of these light rays interacted specularly with the surface according to the approximate Fresnel equations (\ref{eq:Schlicks},\ref{eq:R0}) and the refractive index equations listed in Table~\ref{tab:bestfitequations}. By recording how much of the incident light was reflected we constrained the planar albedo of the surfaces for different compositions and roughnesses. We set the boundary condition that if a light ray lost more than $99.99\%$ of its strength we stopped any further interactions and assumed it was fully absorbed. This was done in order to increase the efficiency of our code.

Fig.~\ref{fig:fractalbouncing} shows an example of a light ray reflecting specularly on our synthetic smoothed, height-corrected magma oceans surface. The bouncing angles may seem distorted because of the scale difference between the x- and y-axes. Fig.~\ref{fig:fractalbouncing} is also a good example of how a light ray may get trapped in a crevice resulting in it bouncing multiple times so a large portion of its original intensity is lost through absorption.
\begin{figure}[h]
	\centering
	\includegraphics[scale=0.80]{Fractal_bouncing.png}
	\caption{An example of a single light ray (red line) interacting with a smoothed height-corrected synthetic magma ocean surface. Take note that the x-axis and y-axis do not have the same scale so the bouncing angle may seem distorted.}
	\centering
	\label{fig:fractalbouncing}
\end{figure}
With the simulator ready we began taking measurements. For each synthetic refractive index ($n_{sy}$) ranging from $1.5-6.0$ we generated multiple fractal terraines for Hurst exponents ranging from $0.95 \: \rm (very \: smooth) \rightarrow 0.01 \: (very \: rough)$. For each Hurst exponent value we simulated 10 morphologically different fractal surfaces where in each case we threw 10,000 light rays and recorded the average reflected component to estimate the planar albedo. However, planets are spherical so the planar albedo could not accurately describe their reflective properties (especially at high latitudes). Because of this, it is necessary to convert the planar albedo into the spherical albedo. By definition, the spherical albedo is the ratio of the monochromatic power reflected by a planet to the monochromatic power incident upon that planet, which is given as follows:
\begin{subequations}
\begin{equation}
    A_{S} = \frac{P_{refl}(\lambda)}{P_{in}(\lambda)}
\end{equation}
Where $P_{refl}$ and $P_{in}$ are the reflected and incident monochromatic powers respectively. The power incident upon the planet is given by:
\begin{equation}
\label{eq:P_in}
    P_{in} = F_{in}(\lambda) \cdot \pi R^{2}_{p}
\end{equation}
Where $F_{in}(\lambda)$ is the monochromatic energy flux and $R_{p}$ is the planetary radius. The power reflected by the planet is more complicated to calculate as one must determine how much light is reflected by each area segment of the planet, as a function of the sub-stellar latitude (or longitude, as for tidally locked planets they are the same). The equation goes as follows:
\begin{equation}
\label{eq:Power_reflected}
    P_{refl} = \int^{\pi/2}_{0} F_{in}(\lambda) \cos{\phi} \times q(\phi) \times 2 \pi R^{2}_{p} \sin{\phi} ~ d\phi
\end{equation}
Where the first section is the monochromatic energy flux received from the star at a given sub-stellar latitude, the second part tells us how much light is reflected at a given sub-stellar latitude, and the third is the area at a given sub-stellar latitude. Although it may seem necessary to multiply Eq.~\ref{eq:Power_reflected} by two in order to account for the whole hemisphere ($\phi = -\pi/2 \rightarrow \pi/2 $), one would also be required to multiply by a half to normalise the integral so that the total power arriving at the hemisphere is $F_{in}(\lambda) \cdot \pi R^{2}_{p}$, hence cancelling each other out. Nevertheless, we can now rearrange as follows:
\begin{equation}
    P_{refl} =  F_{in}(\lambda) \cdot \pi R^{2}_{p} \int^{\pi/2}_{0} 2 \sin{\phi} \cos{\phi} \cdot q(\phi) ~ d\phi
\end{equation}
Using the trigonometric identities,
\begin{equation}
\label{eq:P_refl}
    P_{refl} =  F_{in}(\lambda) \cdot \pi R^{2}_{p} \int^{\pi/2}_{0} \sin{2\phi} \cdot q(\phi) ~ d\phi
\end{equation}
\end{subequations}
We can now divide Eq.~\ref{eq:P_refl} by \ref{eq:P_in} to get the spherical albedo:
\begin{equation}
\label{eq:A_S}
    A_{S} = \int^{\pi/2}_{0} \sin{2\phi} \cdot q(\phi) ~ d\phi
\end{equation}
With the exception of idealised constructs like a Lambertian surface, $q(\phi)$ is never a constant. In order to determine $q(\phi)$, we ran another set of ray-tracing simulations for different refractive indices ($n_{sy}$), Hurst exponents ($H$), and sub-stellar latitudes ($\phi$). The analytic approximation of $q(\phi)$, which incorporates the planar albedo, can be found in Sec.\ref{sec:Results}. By solving Eq.~\ref{eq:A_S} one can determine the spherical albedo, but before this is explored it is necessary to briefly discuss the possibility of atmospheric effects.

\subsection{Atmospheric Effects on the spherical albedo}
\label{sec:atmospheric_effects}

In this paper we focus on how the surface properties of a magma ocean effect the spherical albedo values of airless or near-airless, molten super-Earths. However, one could argue that the vaporised mineral atmospheres resulting from highly irradiated magmas could also affect the spherical albedo values. Fundamentally, this will depend on the composition of the magma and the resultant vaporised atmosphere. For volatile-free magmas, the secondary atmospheres would be thin due to the low refractory nature of rocky species \citep{Schaefer2009,Ito2015}. Furthermore, 1-D atmospheric models show that the high temperatures found on some of these hot super-Earths could prevent the formation of clouds within their daysides \citep{Mahapatra2017}. Therefore, one can consider the case of a thin, cloudless, mineral atmosphere on a hot molten super-Earth. Within this framework, two effects need to be considered. First, the Rayleigh scattering of gases should increase the planetary spherical albedo. Second, the spectroscopic absorption by the gases would decrease it.

The effects of Rayleigh scattering can be understood by considering Eq.~\ref{eq:rs}:
\begin{equation}
\sigma_{\lambda} = \frac{128\pi^5}{3\lambda^4}\alpha^2,
\label{eq:rs}
\end{equation}
where $\sigma_{\lambda}$ is the Rayleigh scattering cross section for a given species, $\lambda$ is the wavelength of light, and $\alpha$ is the static average electric dipole polarisability of the gaseous species being considered. According to the \citet{CRC92}, the $\alpha$ values for the following species (in units of 10$^{-24}$ cm$^3$) are: 24.11 for Na, 43.4 for K, 8.4 for Fe, 5.53 for Si, 10.6 or 11.1 for Mg, 0.802 for O, and 1.5689 for O$_{2}$ \citep[for more information on these species, see][ references therein]{CRC92}. Due to having a large $\alpha$ and being one of the most abundant gases in mineral atmospheres, Na is expected to be the most scattering gas within these extreme environments \citep{Schaefer2009,Ito2015}. The effects of a pure Na mineral atmosphere is illustrated in Fig.~\ref{fig:AS} of Sec.~\ref{sec:Atmosphere} in the appendix where the optical depth from Rayleigh scattering is given for the case of Kepler-10b. It becomes clear how even for a strongly scattering mineral atmosphere, the influence on the spherical albedo becomes negligible for $\lambda\gtrsim$1~$\mu$m.

Conversely, if the magma were to include several volatile elements such as H, C, N, S and Cl, then it would most probably be thick \citep{Schaefer2012}. In that case the scattering of the atmosphere would affect the spherical albedo value. However, the stability of highly volatile elements on close-in exoplanets against photoevaporation is still not well understood \citep{Lopez2017}. For simplicity, in our paper we will only focus on the case where there is a volatile-free mineral atmosphere. However, we are aware that thick, cloud-free, atmospheres would have high albedo values in the optical range \citep{Hamano2015}.

Regarding Kepler-10b, whether the mineral atmosphere affects the albedo or not will depend on its composition, which we do not know with certainty. We give the case of sodium as it is the most optically thick of all non-volatile mineral atmospheric components (and hence, it is an interesting gas to consider). Sodium would affect the albedo at $\lambda= 0.5~\rm \mu m$ as its optical depth is $\simeq 1$. Nevertheless, our code is best adapted for larger wavelengths, primarily $\lambda > 1~\rm \mu m$, so the effects of a mineral atmosphere should not be large. We are aware that the \textit{Kepler} telescope operates at shorter wavelengths so one would expect the mineral atmosphere to contribute towards the spherical albedo. However, the high bolometric Bond albedo values inferred for Kepler-10b cannot fully be accounted for by a mineral atmosphere \citep{Rouan2011}, which is why we argue that an exotic magma must also be present.

We believe that in future works our code could be implemented into radiative transfer models where one could include clouds or a scattering atmosphere (see Sec.~\ref{sec:future_work}). Nonetheless, it is beyond the scope of our study to include atmospheric effects on the planetary spherical albedo as our paper focuses primarily on surfaces.

\section{Results}
\label{sec:Results}

We were able to find an analytic relationship between the planar albedo ($A_{pla}$) of a molten surface and the Hurst exponent (\textit{H}) as shown in Eq.~(\ref{eq:Albedo}) and (\ref{eq:psiconstant}). These equations were found by fitting the numerical data from our simulations. The raw data and its comparison to the planar albedo values predicted by Eq.~(\ref{eq:Albedo}) and (\ref{eq:psiconstant}) can be found in the Supplementary File 1. Within Supplementary File 1, we also show that the $\chi^{2}$ values of our analytic approximation are small and hence we recommend its use over the raw data. For Hurst exponents greater than 0.95 the effects of roughness become negligible so Eq.~\ref{eq:R0} is an appropriate approximation for the planar albedo.
\begin{table}
	\centering
	\caption{Parameters for Eq.~(\ref{eq:Albedo}) and (\ref{eq:psiconstant})}
	\label{tab:parametervalues}
	\begin{tabular}{P{1.5cm}P{1.5cm}}
		\hline
		\hline
		Parameter & $H \leq 0.95$ \\
		\hline
		$\alpha_1$ & -0.397910\\
		$\alpha_2$ & 0.149320\\
		$\alpha_3$ & 0.388928\\
		$\alpha_4$ & 0.852872\\
		$\beta_1$ & 1.594982\\
		$\beta_2$ & -3.375571\\
		$\beta_3$ & 2.437896\\
		$\beta_4$ & -0.881001\\
		$\beta_5$ & 0.194354\\
		$\beta_6$ & -0.041750\\
		$\beta_7$ & -0.081323\\
		\hline
	\end{tabular}
\end{table}
\begin{subequations}
	\begin{equation}
	A_{pla} \approx (\alpha_1 H + \alpha_2 ) + (\alpha_3 H + \alpha_4 )(1 - e^{\Psi \cdot n_{sy}})
	\label{eq:Albedo}
	\end{equation}
	with
	\begin{equation}
	\Psi = \beta_1 H^{6} + \beta_2 H^{5} + \beta_3 H^{4} + \beta_4 H^{3} + \beta_5 H^{2} + \beta_6 H + \beta_7
	\label{eq:psiconstant}
	\end{equation}
\end{subequations}
Where $\alpha_{1\rightarrow4}$ and $\beta_{1\rightarrow7}$ are the coefficients listed in Table~\ref{tab:parametervalues}. Fig.~\ref{fig:Albedo_H_n} is a plot of Eq.~(\ref{eq:Albedo}, \ref{eq:psiconstant}).
\begin{figure}[h]
	\centering
	\includegraphics[scale=1]{Albedo_H_n.png}
	\caption{The planar albedo as a function of the synthetic refractive index ($n_{sy}$) and the Hurst exponent ($H$). This figure was generated with Eq.~(\ref{eq:Albedo}) and (\ref{eq:psiconstant}).}
	\centering
	\label{fig:Albedo_H_n}
\end{figure}
Combining this result with the refractive indices from \citet{Polyanskiy2008} and \citet{Grainger2008} gave us the planar albedo values for different materials at varying roughnesses and wavelengths. To convert the planar albedo into the spherical albedo, the dependence of the albedo on the sub-stellar latitude, $q(\phi)$, was required. From our simulations (see Supplementary File 2 for the raw data and $\chi^{2}$ values) we found that $q(\phi)$ is well approximated as:
\begin{subequations}
\begin{equation}
\label{eq:phase_function}
    q(\phi) \approx 1 - \frac{1 - A_{pla}}{1 + \left(\frac{\phi}{ \mathcal{C}_{1}(H) } \right)^{ \mathcal{C}_{2}(H) }}
\end{equation}
Where,
\begin{equation}
\label{eq:phase_function_constants}
\begin{split}
    &\mathcal{C}_{1}(H) = 0.320717H + 1.14083 \\
    &\mathcal{C}_{2}(H) = 8.54326 + 82.6376H^{3.88648}
\end{split}
\end{equation}
\end{subequations}
From Eq.~\ref{eq:phase_function} it follows that as $\phi \rightarrow \pi/2$, $q(\phi) \rightarrow 1$ because at the terminator the planet is fully reflecting (at right angles there is no absorption). Conversely, when $\phi \rightarrow 0$, the local albedo becomes identical to the planar albedo as the light is interacting with the sub-stellar point. Combining Eq.~\ref{eq:A_S} with Eq.~\ref{eq:Albedo}, (\ref{eq:psiconstant}), (\ref{eq:phase_function}), and (\ref{eq:phase_function_constants}) gives the spherical albedo of a planet. Hence, it is now possible to explore the reflective properties of planets with different magma ocean properties.

We begin by considering the spherical albedo of simple end-member magma compositions (see Fig.~\ref{fig:materials} in the appendix). We adopted a Hurst exponent of 0.8 which is consistent with the type of waves expected on Kepler-10b (see Sec.~\ref{sec:hurst} in the appendix for more information). The best-fit equations for the numerical data of \citet{Polyanskiy2008} and \citet{Grainger2008} are shown in Table~\ref{tab:bestfitequations} of Sec.~\ref{sec:listofequations} in the appendix. A brief summary of our results for the different materials is explained below (for a Hurst exponent of 0.8). The graphs showing the spherical albedo are shown in Sec.~\ref{sec:spectral_albedo_materials} of the appendix.

Planets composed of common Earth minerals and rocks typically have spherical albedo values of $\sim 10\%$. They also tend to drop as the wavelength increases from $\rm 0.5~\mu m$ to $\rm 7.8~\mu m$. The only exceptions are $\rm SiO_{2}$ (i.e. rhyolites) and komatiites which have sudden increases in their spherical albedos at $\rm \sim 7.2~\mu m$.

Pure metal planets would have very high spherical albedos, especially if they were made of Ag, Al, Au, Cu, Fe, Mg, Mo, or Ni as they would reach values close to $100\%$ at larger wavelengths. Metalloids (Ge and Si) would produce lower spherical albedos and within the wavelengths of interest, they appear to stabilise at $30-50\%$.

The reflective properties of metal oxides vary strongly depending on the species. FeO is the shiniest oxidised material analysed in this study with values reaching as high as $45\%$. Conversely, $\rm Al_{2}O_{3}$ and MgO are the least reflective oxidised metallic species with spherical albedo values close to common Earth minerals and rocks such as the ones explained above.

We present two carbonaceous minerals; $\rm CaCO_{3}$ and SiC. SiC's spherical albedo drops from $\sim 20 \%$ to $\sim 15 \%$ from $\rm 0.5-7.8~\mu m$ respectively. Conversely, $\rm CaCO_{3}$ has a spherical albedo $\lesssim 10\%$ for wavelengths below $\rm \sim 6.2~\mu m$ but beyond this it quickly increases to a maximum of almost $\sim 90\%$ at $\rm \sim 6.8~\mu m$. This peak then drops to $\sim 15\%$ for wavelengths greater than $\rm \sim 7.5~\mu m$.

Whilst Fig.~(\ref{fig:materials}) are useful for understanding the spherical albedo spectrum of pure materials; AMOPs are composed of multiple chemical species. It is therefore important to find geochemically self-consistent mixtures of the aforementioned materials so that from an observed albedo spectrum one can categorise the composition of an exoplanet. Hence, we propose the following categorisations for constraining exoplanet compositions from their spherical albedos: evolved molten bulk-silicate Earth (BSE) planets, metallic planets, coreless terrestrial planets (CTPs), reduced planets, and carbon-rich planets.

\subsection{Evolved Magmatic BSE Planet}
\label{sec:evolvedbse}

When a terrestrial planet with a BSE composition is subjected to temperatures above $\rm\sim1500~K$, the minerals found within the magma begin to vaporise \citep[e.g.][]{Schaefer2009,Schaefer2012,Ito2015}. This vaporisation could lead to a strong evolution in the chemical composition of the magma \citep{Kite2016}. It therefore follows that the physical and optical properties of the lava should also change with time. In Fig.~\ref{fig:evolvedterrestrial} we show our predicted spherical albedo spectrums for a BSE airless planet under 4 different levels of vaporisation: $0\%$, $\lesssim80\%$, $\sim 92\%$, and $\sim 100\%$. The first has a typical BSE composition \citep{Oneill1998} whilst the latter three approximately correspond to $60:40$ of $\rm MgO-SiO_{2}$, $50:50$ of $\rm CaO-Al_{2}O_{3}$, and pure $\rm Al_{2}O_{3}$ respectively.

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.9]{evolved_BSE.png}
	\caption{The spherical albedo for a terrestrial planet with an initial composition corresponding to the BSE that has been evolved due to chemical interactions between the magma and the mineral atmosphere. Four different spherical albedo spectrums are shown corresponding to varying levels of vaporisation. The BSE composition is from \citet{Oneill1998} whilst the vaporised compositions were taken from \citet{Kite2016}. The light-blue section is \textit{Kepler's} band-pass. This figure was generated with Eq.~(\ref{eq:A_S}), (\ref{eq:Albedo}), (\ref{eq:psiconstant}), (\ref{eq:phase_function}), and (\ref{eq:phase_function_constants}).}
	\centering
	\label{fig:evolvedterrestrial}
\end{figure}

Our synthetic spectrums predict spherical albedos lower than $\sim 11\%$. Due to current limitations in our astronomical instrumentation there appears to be no clear way of distinguishing between different types of Earth-like compositions regardless of how evolved the magma is. Notwithstanding, our findings suggest than a spherical albedo value of $\sim 11\%$ or lower, with a decrease at longer wavelengths, is suggestive of a composition similar to the Earth's.

\subsection{Metallic Planet}
\label{sec:metallic}

Prior to considering the spherical albedo values of metallic planets it is important to understand their formations. One can begin by analysing Mercury that has an estimated bulk iron abundance of $\sim70\%$ \citep[e.g.][]{Riner2008,Malavergne2010,Hauck2013,Chabot2014}. There are three main explanations for this super-ferruginous composition: thermally-induced surface vaporisation \citep[e.g.][]{Cameron1985}, collisional stripping \citep[e.g.][]{Benz1988}, and chemical/thermal gradients within the primordial solar nebula triggering an unusual formation \citep[e.g.][]{Lewis1972,Lewis1974}. For larger bodies such as super-Earths, surface vaporisation would be unable to trigger sufficient mass-losses as the required escape velocities of silicate species are too large (Ito et al. in prep). Regarding collisional stripping, \citet{Marcus2010} showed that the maximum amount of iron that a super-Earth could have is $\sim70\%$ which is compatible with the abundance in Mercury. Finally, forming iron-rich exoplanets due to chemical/thermal gradients has been shown to be possible in our solar system \citep[e.g.][]{Lewis1972,Lewis1974}, and considering how iron is the sixth most abundant element in the universe (by mass) there does not appear to be a strong argument against this mechanism either. We are therefore left with the latter two mechanisms for the formation of metal-rich super-Earths.

Before one can explore the outcomes that could arise from these two dissimilar formation mechanisms, it is necessary to issue the following caveat: a bulk iron composition is not implicative of a surface also containing iron. For example, Mercury's surface composition is very iron-poor $\rm 1-2~wt.\%$ \citep{Evans2012,Weider2014} which may seem paradoxical at first. There is a straightforward explanation for this; Mercury formed from very reduced materials. We will discuss reduced planets more in-depth in section \ref{sec:reduced}. Notwithstanding, if a terrestrial exoplanet were to form under very oxidised conditions, a large portion of its total iron would bind with oxygen to form the usual oxides (FeO, $\rm Fe_{2}O_{3}$, and $\rm Fe_{3}O_{4}$). These compounds are less dense than pure metals or alloys so they would become incorporated into rock-forming minerals such as olivine and pyroxenes that would then form the mantle and crust \citep{Elkins2008}. Conversely, the pure iron and iron-nickel alloys would sink to the centre of the planet thus forming a core. Following this line of thought, after a certain critical level of oxygen, all of the iron would form oxides therefore leaving the planet coreless (see Sec.~\ref{sec:coreless}). It is therefore very important not to assume that an iron-rich world will host an iron-rich surface as the stability of metals in the crust and mantle are strongly dependent on the presence of oxygen. Keeping in mind the above caveat, one can now refer back to the two distinct formation mechanisms and show how they can lead to two geochemically distinct magmas.

According to the collisional-stripping model, a planet becomes Fe-rich \textit{after} its formation due to the partial removal of the mantle, while its core remains relatively intact. Therefore, the abundance of iron in the mantle and crust (from which the magmas are formed) is approximately independent of any collisions (excluding exogenous deposits). It then follows that the Fe-richness of metallic exoplanets' magmas is more strongly dependent on the presence of oxygen than the bulk abundance of iron. Conversely, if iron planets form due to thermal/chemical gradients in the primordial stellar nebula then there could be a relationship between being oxygen-poor and iron-rich. This is best explained by briefly summarising the requisites for both of these outcomes:
\begin{itemize}
	\item At very high temperatures $\rm MgSiO_{3}$ is weakly captured whilst iron condenses with relative ease \citep[e.g.][]{Lewis1972,Lewis1974}. Consequently, under these conditions the minerals that form are richer in iron-bearing minerals and poorer in silicates than elsewhere in the stellar nebula.
	\item At very high temperatures the chondrules that form enstatite chondrites are created from supercooled, fluid, silicate, rain-like droplets \citep{Blander2009}. These reduced chondrites are believed to be the materials that form reduced planets such as Mercury \citep[e.g.][]{Nittler2011,Malavergne2014,Nittler2019}.
\end{itemize}
In contrast to the collisional stripping model, if iron-planets form from chemical/thermal gradients, having a high abundance of iron might coincide with being reduced. This would result in most of the iron being captured within the core and little being left within the mantle. Therefore, according to this mechanism, metallic exoplanets could have felsic magmas that are sulfur-rich \citep[e.g.][]{Haughton1974,Mccoy1999,Namur2016}.

To summarise, the magma composition of metallic planets is strongly dependent on the formation mechanisms and history of the planet. Hence, we believe that the bulk iron mass-fraction cannot be constrained from the measured spherical albedo spectrum of these types of bodies.

\subsection{Coreless terrestrial Planet}
\label{sec:coreless}

In order to model coreless terrestrial planets (CTPs) we first had to set a planetary bulk composition and then subject it to the relevant thermodynamic conditions. In order to do so we used the elemental composition listed in Table~\ref{tab:chemicalcomposition} with the temperature and pressure shown in Table~\ref{tab:Hurstcalculation} (see appendix) and then applied the GG\textsc{chem} geochemical code \citep{Woitke2018}. Our resultant magma composition is well approximated as $60:40$ of iron oxides to enstatite which is much more iron-rich than the BSE ($\rm \lesssim 10\%$ iron oxides). However, because GG\textsc{chem} simulates only FeO as the Fe-bearing oxide in liquid phases, it does not specify the ratio of $\rm FeO:Fe_{2}O_{3}:Fe_{3}O_{4}$. As shown in Fig.~\ref{fig:materials}, the spherical albedo values of these species differ moderately so one cannot treat all of them as FeO. We will therefore present the upper- and lower-bounds with the average spherical albedo for a given Hurst exponent (Fig.~\ref{fig:corelessterrestrial}).

\begin{table}
	\centering
	\caption{Predicted Bulk Chemical Composition of Coreless Terrestrial Planets (CTPs)}
	\label{tab:chemicalcomposition}
	\begin{tabular}{P{2cm}P{3cm}}
		\hline
		\hline
		Chemical & Abundance ($\%$) \\
		\hline
		O &	35.2\\
		Si & 13.5 \\
		Fe & 37.8\\
		Mg & 11.3\\
		Ca & 1.1\\
		Al & 1.1 \\
		\hline
	\end{tabular}
	\tablecomments{Based on the findings from \citet{Elkins2008}}
\end{table}

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.9]{coreless_terrestrial.png}
	\caption{The spherical albedo of a coreless terrestrial planet with the composition given by \citet{Elkins2008}. The red dotted line is the average coreless spherical albedo value for $H = 0.8$, the red shaded region represents the uncertainty, and the black line is for the spherical albedo value of a cored BSE planet \citep{Oneill1998}. The light-blue section is \textit{Kepler's} band-pass. This figure was generated with Eq.~(\ref{eq:A_S}), (\ref{eq:Albedo}), (\ref{eq:psiconstant}), (\ref{eq:phase_function}), and (\ref{eq:phase_function_constants}).}
	\centering
	\label{fig:corelessterrestrial}
\end{figure}

Our results show that coreless bodies have much higher spherical albedos than their differentiated counterparts with similar compositions. This is due to iron-oxides being shinier (see Fig.~\ref{fig:materials}) and more abundant in CTP magmas. The uncertainty in the spherical albedo is due to the different possible ratios of $\rm FeO:Fe_{2}O_{3}:Fe_{3}O_{4}$, which is dependent not only on the Gibbs energies but also on how much oxygen was present during the formation of the body. Whilst one might not be able to differentiate between which iron oxides are present within the magma, one could make a strong argument for the presence of iron within these bodies. Adjusting the surface roughness parameter (\textit{H}) and maximising the FeO abundance allows for a theoretically maximum spherical albedo of $\sim 55\%$.

\subsection{Reduced Planet}
\label{sec:reduced}

In order to understand the formation mechanism of reduced planets, one can analyse Mercury which is the most reduced planet in our solar system \citep[e.g.][]{Nittler2011,Malavergne2014,Nittler2019}. This is probably because it formed from enstatite chondrites which are extremely reduced \citep{Blander2009}. These rocks can be classified into two groups; EH (high enstatite) and EL (low enstatite) chondrites. Each group has slightly different chemical makeups \citep{Wasson1988} but their close similarities allow for an approximate average composition to be deduced \citep[see Table~\ref{tab:chemicalcomposition2},][]{Javoy2010}.
\begin{table}
	\centering
	\caption{Predicted Bulk Chemical Composition of Reduced Planets}
	\label{tab:chemicalcomposition2}
	\begin{tabular}{P{2cm}P{3cm}}
		\hline
		\hline
		Chemical & Abundance ($\%$) \\
		\hline
		O &	33.3 \\
		Si & 20.2 \\
		Fe & 28.3 \\
		Mg & 14.0 \\
		Ca & 1.0 \\
		Al & 0.97 \\
		Ni & 1.73  \\
		Ti & 0.06 \\
		Cr & 0.36 \\
		Co & 0.08 \\
		\hline
	\end{tabular}
	\tablecomments{Based on the findings from \citet{Javoy2010}}
\end{table}

By applying the GG\textsc{chem} code \citep{Woitke2018} with the elements listed in Table~\ref{tab:chemicalcomposition2}, the temperature and pressure shown in Table~\ref{tab:Hurstcalculation}, we predict a magma composed of over $\gtrsim 90\%$ pure Fe with the rest mostly comprising typical terrestrial silicate minerals such as forsterite and spinel. The reason why the Fe concentration in the magma is much higher than the one shown in Table~\ref{tab:chemicalcomposition2} is because one also has to consider the composition of the emitted gas. This vaporised atmosphere contains many of the volatile elements that used to be present in the melt. Nevetheless, because of gravity, the pure iron would sink to the planet's core thus resulting in an iron-poor mantle. These predictions are strongly corroborated by observational \citep[e.g.][]{Evans2012,Weider2014} and geophysical \citep[e.g.][]{Riner2008,Hauck2013} data of Mercury. Furthermore, spinel, fosterite and many of the other silicate minerals that are predicted by the GG\textsc{chem} geochemical code are the same materials that makeup large portions of Earth's mantle. These silicates would be molten within the high temperature ranges we are focusing on in this paper. In other words, from a reduced planet made of enstatite chondrites (or from similarly reduced rocks) we would expect silicate-rich magmas that are iron-poor. A rhyolitic magma could be a good contender due to being very $\rm SiO_{2}$-rich \citep[typically $\gtrsim 69\%$;][]{Lemaitre2002}. If rhyolites were to be the most common magmas on reduced super-Earths then from Fig.~\ref{fig:materials} it can be seen that spherical albedo values less than $10\%$ would be expected. However, other exotic magmas that are oxygen-poor could also be possible. For instance, if a planet were carbon-rich then carbonaceous magmas might exist. We discuss this concept in Sec.~\ref{sec:carbonplanet}.

Finally, we would like to briefly discuss one potential way of forming a reduced, airless, molten planet with pure iron deposits on its surface. Suppose a planet were to form under very cold conditions, this would result in most of the oxygen bonding with silicon, carbon, hydrogen and magnesium. Considering how iron is not a very aerophilic element, some of it may be left in its elemental form. If this planet were then to migrate very close to its host star, XUV-irradiation may fully erode away its atmosphere leaving behind a denuded core \citep[e.g.][]{Lecavelier2007, Ehrenreich2011, Lammer2013, Owen2013, Jin2014, Owen2017, Jin2018, Kubyshkina2018(2),Locci2019,Modirrousta2020b}. Due to iron being relatively heavy, it may be resilient to hydrodynamic escape \citep[e.g.][]{Zahnle1986,Hunten1987,Luger2015} which could result in pure iron deposits being left behind on the bare core. Pure iron is very reflective (see Fig.~\ref{fig:materials}) so it is possible for these planets to have very high spherical albedo values if there were to be a copious abundance of iron in its elemental form. Undeniably, the stability of pure iron within a lava matrix would be subjected to the usual caveats of the circulation and sinking timescales. These systems are highly speculative and warrant more research.

In summary, the spherical albedos of reduced planets that form from enstatite chondrites should be very low ($\lesssim 10\%$) as their surfaces are expected to be rich in silicate species. We find that the amount of light they reflect is similar to that of BSE composition planets, albeit with slighly different spectrum shapes. However, one cannot rule out the possibility of more unusual formation mechanisms which could lead to an ambiguity in the predicted spherical albedos.

\subsection{Carbon-rich Planet}
\label{sec:carbonplanet}

There are several arguments for carbon-rich super-Earths such as the existence of carbonaceous bodies in our solar system (e.g. c-type asteroids) and the presence of main sequence stars with high carbon-to-oxygen ratios (C/O). These stars could host these exotic worlds, especially since theoretical modelling predicts that even with moderate C/O ratios several carbon-rich worlds could form \citep{Moriarty2014}. However, as of February 2020 there are no known carbon-rich super-Earths. Because of this limitation we will base our magma composition on theoretical principles. We will avoid using pure carbon \citep{Madhusudhan2012} as the minimum melting points of carbon allotropes tend to be $\rm > 4000~K$ \citep[e.g.][]{Bundy1989,Bundy1996,Correa2006}, which is generally beyond the equilibrium temperature of most super-Earths. Hence, we will adopt a more chemically diverse composition such as the one delineated in \citet{Miozzi2018}.

Whilst SiC has been found in presolar grains \citep[e.g.][]{Lodders1995,Hoppe2010}, it has a tendency to decompose due to the high affinity between oxygen and silicon. For instance, if a planet were to be composed of equal parts Si-O-C, most of the silicon would bind with the oxygen to form silicates such as $\rm SiO_{2} $ and SiO. The carbon would then bond with itself to form graphite in the outer layers of the planet and diamond at greater depths. Due to graphite being less dense than $\rm SiO_{2} $ and having a higher melting temperature, it would rise to the top layers of the mantle in the form of solid graphitic grains. Our simple analogy is strongly supported by high-pressure, high-temperature (high P-T) laboratory experiments \citep[e.g.][]{Hakim2018,Hakim2019} that arrived at the same conclusion, albeit in more detail. These high P-T experiments predict that extremely reduced conditions would be required in order to maintain SiC in a geochemically stable condition. Therefore, in spite of SiC planets being theoretically plausible, they may be unlikely. Notwithstanding, we included the spherical albedo of SiC in Fig.~\ref{fig:materials}.

On that account, it seems more probable that carbon-rich worlds will be composed of silicate mantles covered in solid graphite. The spherical albedo values of these planets would therefore be dominated by the graphitic layer which we cannot model with our code due to it being adapted for fluid systems. If one ignores surface roughness, the spherical albedo of a graphite layer within the range of $\rm 0.5-7.8~\mu m$ would vary from $30-80\%$ \citep[e.g.][]{Taft1965,Philipp1977,Querry1985,Djurisic1999,Kuzmenko2008,Papoular2014}, which would make carbon-rich super-Earths the brightest airless bodies listed in our paper.

\section{Discussion}
\label{sec:Discussion}

Measuring, constraining, and interpreting the spherical albedo of exoplanets is an active area of research which is mostly limited by observational constraints. With the JWST expected to launch in 2021, the \textit{Atmospheric Remote-sensing Infrared Exoplanet Large-survey} (ARIEL) in 2028, a more in-depth analysis could be possible. Because of this we now provide two tests that could be used to constrain the magma composition of molten super-Earths from their spherical albedo values:

\subsection{Optimistic Test: High Precision High Accuracy Measurements (HPHA)}

In the ideal case scenario, the extracted albedo spectra from the aforementioned missions would be robust enough to perform statistical tests on them (e.g. MCMC algorithms). The shape of the spectrum would be the major source of information for the surface composition whilst the amplitude would mostly determine the surface roughness (i.e. the Hurst exponent). The experimental data that relates a material's refractive index to the wavelength of light began almost a century ago and it is ongoing. From these relationships, Eq.~(\ref{eq:nsynth}), (\ref{eq:A_S}), (\ref{eq:Albedo}), (\ref{eq:psiconstant}), (\ref{eq:phase_function}), and (\ref{eq:phase_function_constants}) could be used to constrain the surface compositions of molten, airless super-Earths. Because of this, regardless of whether ARIEL or the JWST missions could extract HPHA measurements, it is certainly possible that in the future the spherical albedo spectrum of a planet could be matched to a certain composition in a similar fashion to how exoplanet atmospheres are being constrained from spectroscopic data.

\subsection{Conservative Test: Low Precision Low Accuracy Measurements (LPLA)}

A more realistic scenario is that a few albedo measurements at varying wavelengths with considerable uncertainties are extracted. For this scenario we propose a test that will most probably be possible for ARIEL and the JWST missions once an optically thick atmosphere is excluded (i.e. the planet is airless or near-airless). We propose that based on the mass and radius of an observed AMOP one limits the possible internal compositions. From this information, the observed albedo spectrum is compared to our five compositional categories described in sections \ref{sec:evolvedbse}, \ref{sec:metallic}, \ref{sec:coreless}, \ref{sec:reduced}, and \ref{sec:carbonplanet}. For instance, in the case of Kepler-10b the mass and radius measurements are consistent with a coreless or near-coreless interior which is illustrated in Fig.~\ref{fig:kepler10b}. Therefore, whilst the spherical albedo is very uncertain, an extremely ultramafic magma is compatible with the values shown in Fig.~\ref{fig:corelessterrestrial}.

\subsection{The Effects of Changing the Hurst Exponent (H)}

In Fig.~\ref{fig:Albedo_H_n} we show how the Hurst exponent and the synthetic refractive index affect the planar albedo of a system. However, understanding the isolated influence that the \textit{H} parameter has may be difficult to picture so we have included Fig.~\ref{fig:FeOH}. Although the shape of the spectrum is slightly affected by varying \textit{H}, the most prominent change is in the amplitude. This is an important point to consider as it implies that strong deviations in the shape of the retrieved albedo spectrum are most probably due to compositional differences between the adopted (predicted) magmas and the ones that are observed. From an observational perspective, this lowers the degeneracy of the inferred surface properties from the observed spherical albedo as the effects from the surface roughness (\textit{H}) and the composition can be partially isolated.
\begin{figure}[h]
	\centering
	\includegraphics[scale=0.8]{Albedo_FeO_H.png}
	\caption{The spherical albedo of FeO for different Hurst exponents and wavelengths. The light-blue section is \textit{Kepler's} band-pass. This figure was generated with Eq.~(\ref{eq:A_S}), (\ref{eq:Albedo}), (\ref{eq:psiconstant}), (\ref{eq:phase_function}), and (\ref{eq:phase_function_constants}).}
	\centering
	\label{fig:FeOH}
\end{figure}

\subsection{Comparing our Results with \citet{Hu2012}}

Our study focuses on molten airless bodies whilst \citet{Hu2012} focuses on solid ($\rm T_{eff} = 300 - 880$ K) planetary surfaces so there will be inevitable differences in our results. For instance, for solid systems the surface roughness is affected by grain configurations and sizes \citep[e.g.][]{Bennett1961}, which is very hard to model mathematically. In their study, fine grained surfaces were assumed when constructing their synthetic spectra. This will lead to an overestimation of the albedo value as finer grains typically have higher albedos than their coarser-grained counterparts. Because of this they point out that their results are “upper limits”. This is further emphasised in our study, as in Fig.~\ref{fig:FeOH} we show that differences in the surface roughness could change the spherical albedo values by $\sim 15\%$ for certain materials. However, this alone might not fully explain the different albedo values calculated by our study and theirs, so below we will explore other possibly reasons:
\begin{itemize}
	\item\textbf{Basaltic:} Their basaltic values are much higher than our results as well as \citet{Pollack1973}, \citet{Egan1975}, \citet{Shestopalov2013}, and \citet{Fornasier2016}. One possibility is that they used very young basaltic powder which can be several times more reflective than older specimens. The reason why there is a darkening with age is due to interactions between the basalt grains and the H-ions from solar winds \citep{Shestopalov2013,Fornasier2016}.
	\item\textbf{Fe-oxydised}: In their paper, they adopt a composition of $50\%$ $\rm Fe_{2}O_{3}$ (haematite) and $50\%$ basalt, with which they predict an albedo of $10-20\%$. This is consistent with our results as we predict a basaltic albedo of $2-7\%$ and a value of $20-30\%$ for $\rm Fe_{2}O_{3}$, so a $50:50$ mixture could reproduce their results.
	\item\textbf{Metal-rich:} We consider planets that are rich in metals (bulk), which as explained in Sec.~\ref{sec:metallic} would result in metal-poor surfaces. Conversely, in \citet{Hu2012} they adopt an $\rm FeS_{2}$ composition.
	\item\textbf{Ultramafic:} Their ultramafic albedo is large, but it could be explained by a very high concentration of metal oxides (see Sec.~\ref{sec:coreless}). However, if we assume that \citet{Hu2012} used a moderate metal oxide abundance then their results would be substantially larger than ours (see Fig.~\ref{fig:corelessterrestrial}). This may be explained by their fine-grained assumption and the possibility that they used very young basaltic specimens.
\end{itemize}
Whilst the results from \citet{Hu2012} differ moderately from ours, we are working in two different regimes: molten versus solid. Thus, we do not feel like there is a contradiction between our spherical albedo values and theirs. However, we are aware that our model has its own limitations that could be improved upon in the future. We list these limitations in Sec.~\ref{sec:limitations}.

\subsection{Limitations of our Model}
\label{sec:limitations}

Our model has limitations that need to be noted by anyone who wishes to use it. Here we present a list of shortcomings and how they could be improved upon by future research:
\begin{itemize}

    \item\textit{We used Earth's oceans as an analogue for the roughness of the magma oceans of airless or near-airless molten planets.} We used this approximation twice: when applying the Gaussian filter to the pre-smoothed synthetic magma ocean surfaces (i.e. Fig.~\ref{fig:fractalsurfacesharp}), and when solving for the constant $h_{0}$ in Eq.~\ref{eq:DurstH}. The question one might ask is whether the viscosity of lava is comparable to that of water in order for surface perturbations (such as waves) to be similar. Using the Vogel–Fulcher–Tammann equation \citep{Vogel1921,Fulcher1925,Fulcher1925b,Tammann1926} the similar dynamical viscosities of these two fluids can be shown:
    \begin{equation}
        \log_{10}\left(\mu \right) = A + \frac{B}{T - C}
    \end{equation}
    Where $\mu$ is the dynamic viscosity at $T$ which is the temperature of interest; $A$, $B$, and $C$ are experimentally-determined constants, which depend on the composition of the magma. If one were to set the temperature high enough one can see how the viscosity drops precipitously. Let $T$ be the day-side temperature of Kepler-10b $\approx 2750~\rm K$ and for the experimental constants we will adopt the values for peridotite \citep{Dingwell2004}: $A \simeq -4.31$, $B \simeq 3703$, $C \simeq 761.7$. These values give $\mu \approx 10^{-2.5}~\rm Pa~s$ which is only twice as viscous as water \citep{Viswanath1989}. This result is compatible with other studies which are consistent with water-like lava viscosities at very high temperatures \citep[e.g.][]{Russell2004,Giordano2008}. A similar calculation can be done for metals to show that at high temperatures they have very low viscosities \citep[e.g.][]{Alfe1998}.

    It is important to note that our model will only work for fluid surfaces. For example, for silicate planets effective temperatures $\gtrsim 1500~\rm K$ will result in magma oceans. Moreover, higher temperatures will increase the accuracy of our model as the melt viscosities will approach that of water. For other materials the conditions will differ, such as a carbonatite magma that at $\simeq 800~\rm K$ will be molten and have a viscosity similar to that of water. Nevertheless, we are aware that the viscosity of the magma oceans depends on many parameters such as the composition, the spherical albedo, and the effective temperature of the planet. These are important parameters that should be investigated further but they are beyond the aims of our study.

	\item\textit{We assumed that the refractive index is independent of the temperature of the system}. We are aware that this is a strong simplification as, for example, Ge and Si's refractive indices grow by approximately $\sim 5\%$ as their temperatures increase from 100 K to 600 K \citep{Li1980}. We also lack experimental data on how the refractive indices evolve when there is a phase change which could be problematic as we use the refractive indices for the solid phases whilst magmas are fluids. However, as shown by \citet{Essack2020}, the albedo of $100\%$ molten and $100\%$ solidified basaltic glass remains relatively constant which implies that the refractive indices might not vary too strongly. In any case, the dependence of the refractive indices on the temperature and phase of the materials are poorly constrained and warrant more research. This limitation can be fixed with relative ease; whenever a better data set is discovered for a given material one can replace the $n$ and $k$ or $n_{sy}$ values given in table~\ref{tab:bestfitequations} with the improved ones. This new data could then be used within Eq.~(\ref{eq:A_S}), (\ref{eq:Albedo}), (\ref{eq:psiconstant}), (\ref{eq:phase_function}), and (\ref{eq:phase_function_constants}).

	\item\textit{We did not account for the albedo contribution from small, solid, undissolved particles i.e. single-scattering albedo.} The presence of solid undissolved particles could strongly influence the spherical albedo of a planet. As explained by \citet{Rouan2011}, molecules such as $\rm UO_{2} $ and $\rm ThO_{2} $ would not dissolve in certain melts such as an $\rm Al_{2}O_{3}-CaO $ matrix and they would be solid. Other molecules such as MgO and CaO also have very high fusion temperatures but they are more commonly dissolved. Notwithstanding, spherical albedos as high as $\sim 50\%$ may be caused by such minerals, so planets like Kepler-10b may be very reflective because of their presence. Our model does not account for these particles but one could separately include their contribution therefore overcoming this limitation.

	\item\textit{We assume that the spherical albedo becomes independent of the Hurst exponent for $\rm H>0.95$.} In reality the spherical albedo would change for surface roughnesses ranging from 0.95 to 1.0. However, it becomes exceedingly difficult to model such surfaces as the ocean stops behaving fully like a fractal and starts having properties more akin to simple geometrical shapes. Because of this, the Davies-Harte \citep{Davies1987} method cannot be used. Hence, considering how the spherical albedo values of surfaces with Hurst exponents of 0.95 are very close to that of $H = 1.0$; we approximated the surfaces as being flat for $0.95 < H < 1.0$ (i.e. Eq.~\ref{eq:R0}).

\end{itemize}

\subsection{Future Work}
\label{sec:future_work}

There is a bimodal distribution of exoplanet radii with one peak at $\rm \sim 1.3R_{\oplus}$, the other at $\rm \sim 1.75R_{\oplus}$, and the minimum at $\rm \sim 2.4R_{\oplus}$ \citep{Fulton2017}. According to theoretical modelling \citep[e.g.][]{Owen2013,Owen2017,Jin2018,Modirrousta2020b} and observational data \citep[e.g.][]{Swain2019} most planets in the first peak are lacking a hydrogen atmosphere whilst most planets in the second peak have large hydrogen envelopes. On the one hand, for planets with hydrogen atmospheres it is much easier to constrain their compositions as atmospheric species can be detected \citep{Tsiaras2016,Ridden2016,Esteves2017} and then used to infer some interior properties. On the other hand, if a planet lacks a large hydrogen component then spectroscopy may not be feasible which leaves one with less data to deduce interior properties. In the case of molten, airless super-Earths (located mostly within the first peak) the spherical albedo could be used in an analogous manner to spectroscopic data in order to deduce the surface composition and hence constrain the interior structure. Hence, one could test exotic bodies like 55 Cancri e which could either be a coreless body \citep{Bourrier2018(2)}, a carbonaceous planet \citep{Madhusudhan2012,Miozzi2018}, an icy planet \citep{Zeng2013,Zeng2016}, or a planet with a hydrogen atmosphere partially confined to the nightside and with a dayside that is near-airless \citep{Modirrousta2020}.

In addition, due to our code being analytic, it could be implemented into atmospheric general circulation models (GCMs) or spectral retrieval codes. This might provide a more holistic analysis of a planet since the effects of an atmosphere and those of the surface could be better coupled. Furthermore, it would be interesting to expand our work to non-molten surfaces such as desert planets. In theory, solid surfaces could also be explained using fractal mathematics \citep[e.g.][]{Mandelbrot1975,Shelberg1983}. We would also welcome more experimental data on the refractive indices of materials at different temperatures, pressures, and phases in order to better model the spherical albedo values of exoplanetary surfaces.

\section{Conclusions}
\label{sec:conclusion}

In this paper we present an analytic model for the spherical albedo of molten, airless or near-airless super-Earths. We developed this model by fitting the results from the simulations of photons interacting with magma ocean surfaces that had varying reliefs and geochemical compositions. Our list shown in Table~\ref{tab:bestfitequations} is not comprehensive as there are many potential chemicals that could be added to the collection. Due to being analytic in nature, our model could be adopted in a computationally inexpensive manner in order to constrain the surface composition of some exoplanets.

\section*{Acknowledgements}

We acknowledge the support of the ARIEL ASI-INAF agreement n. 2018-22-HH.0. We also thank the researchers and staff at the Telescopio Nazionale Galileo, La Palma (Spain) for introducing us to Kepler-10b that inspired us to carry out this study. We are grateful to Podolak M. and the anonymous referees for their useful comments.


\section{Introduction}

The spectrum of the cosmic microwave background (CMB) carries unique
information/constraints on the thermal history of the universe since energy
releases occurring over many redshift decades can leave their imprint on it
\citep{ZeldovichSunyaev1969, SunyaevZeldovich1970, IllarionovSunyaev1975,
DaneseDeZotti1977, Burigana1991b, Chluba2012, KhatriSunyaev2012,
ChlubaJeong2014, Tashiro2014, Chluba2016, DeZotti2016, Chluba2019}. However, at
very high redshifts such imprints are erased by thermalization effects due to
the combined action of photon emission processes and of Compton scattering.
Small distortions are completely thermalized at $z > \hbox{few}\times 10^6$
\citep{DaneseDeZotti1982, Burigana1991b, HuSilk1993a, Chluba2014a}. The
thermalization is less efficient for large distortions which can keep some
visibility up to redshifts higher by factors of several \citep{Chluba2020}.

At still higher redshifts, releases of very large amount of energy could have
occurred without leaving any visible track in the CMB spectrum. On the other
hand, the outcomes of primordial nucleosynthesis (or big-bang nucleosynthesis,
BBN) are sensitive to the radiation energy density. Indeed, just the
consideration of the production of light elements in the early universe led to
the prediction of the CMB \citep{Gamow1948, AlpherHerman1949}.

The present-day accurate determinations of cosmological parameters entail
strong constraints on the CMB energy density both at the BBN epoch and at
electron-photon decoupling, hence on the additional amount of energy that could
have been released after the BBN, at high enough redshifts to ensure a tight
coupling between electrons and radiation. Such constraints are quantified in
Sect.~\ref{sect:results}. Our conclusions are summarized in
Sect.~\ref{sect:conclusions}.

\begin{figure}

	\includegraphics[width=\columnwidth]{CMB_constraints.png}
    \caption{Upper limit (95\% confidence level) to the fractional amount of energy
    that could have been added to the CMB after the primordial nucleosynthesis epoch, derived in this paper (black horizontal line),
    compared with constraints derived by \citet{Chluba2020}. The red dotted horizontal lines show the limits
    obtained from the effective number of neutrino species, assuming that a
fraction $f_\nu$ of the total injected energy goes into neutrinos. The lower and upper lines refer to
$f_\nu=0$ and $f_\nu=0.5$, respectively; note that the value for $f_\nu=0$ has been corrected
as mentioned in the text. Constraints
    from CMB spectral measurements, as computed by \citet{Chluba2020}, are also shown.
    The upper (blue) and lower (pink) dot-dashed lines show the 95\% confidence upper limits
    implied by CMB spectral measurements from COBE/FIRAS ($\delta\epsilon/\epsilon_\gamma< 6\times 10^{-5}$)
    and from a future PIXIE-like experiment ($\delta\epsilon/\epsilon_\gamma< 10^{-8}$),
    respectively. The dotted lines deviating from the dot-dashed lines at their right end, again from \citet{Chluba2020},
    show the corresponding limits obtained under the small-distortion approximation.}
    \label{fig:constraints}
\end{figure}

\section{Constraints on the radiation energy density at
BBN}\label{sect:results}

The results of BBN calculations can be presented as a function of the
present-day dimensionless baryon-to-photon number density ratio $\eta =
n_b/n_\gamma$ \citep{Steigman2007, Cyburt2016, Fields2020}. The BBN redshift
($z_{\rm BBN}\sim 3\times 10^8$) is much higher than the CMB thermalization
redshift even in the case of strong deviations from equilibrium
\citep{Chluba2020}. Hence the CMB photon number density is $n_{\gamma,\rm BBN}
=20.28\, T_{\rm BBN}^3$ and the energy density is $\epsilon_{\gamma,\rm BBN} =
a\, T_{\rm BBN}^4$, $a=7.5657\times
10^{-15}\,\hbox{erg}\,\hbox{cm}^{-3}\,\hbox{K}^{-4}$ being the black-body
radiation density constant and $T_{\rm BBN}$ the CMB temperature at BBN.

The present-day baryon number density $n_{b,0}$ is related to the baryon
density parameter $\omega_b=h^2 \Omega_b$ by;
\begin{equation}
\omega_b = \frac{\langle m_{\rm b}\rangle \ n_{\rm b,0}}{h^{-2} \rho_{\rm crit,0}} \\
\end{equation}
where $h=H_0/100\,\hbox{km}\,\hbox{s}^{-1}\,\hbox{Mpc}^{-1}$ and $\Omega_b=
\rho_b/\rho_{\rm crit} = 8\pi G \rho_b/3 H_0^2$, $G$ being the gravitational
constant and $\rho_b = \langle m_b\rangle n_b$. The mean mass per baryon is
very well approximated by \citep{Fields2020}  $\langle m_b\rangle = (1+\delta)
m_{\rm H}$ where $\delta=-\left[ 1.744  +  7.119  \left( Y - 0.245 \right)
\right] \times 10^{-3}$ accounts for the reduction of the mass due to helium
binding; $Y$ is the primordial helium mass fraction. Combining the above
relations we get \citep{Fields2020}:
\begin{equation}\label{eq:eta}
\eta_{\rm BBN} = \frac{273.754\times 10^{-10}}{1 - 7.131 \times 10^{-3} \left(Y-0.245
\right)} \, \left(\frac{2.7255\rm K}{T_{\rm BBN,0}}\right)^3 \,\omega_b\, ,
\end{equation}
where $T_{\rm BBN,0}$ is $T_{\rm BBN}$ redshifted to $z=0$. Note that this
expression does not contain the Hubble parameter and therefore is not affected
by the current discrepancy between the value derived by \textit{Planck} CMB
anisotropy measurements \citep{Planck2018parameters} and the local value
derived by \citet{Riess2019}.

The excellent agreement between the values of $\eta$ inferred from BBN and from
CMB anisotropies \citep{Planck2018parameters, Fields2020} implies a tight upper
limit to the amount of energy that could have been released to the CMB after
the nucleosynthesis epoch. \citet{Fields2020} found $\eta_{\rm BBN}=(6.084\pm
0.230)\times 10^{-10}$ and $\eta_{\rm CMB}=(6.090\pm 0.060)\times 10^{-10}$,
allowing for variations of the number of neutrino species ($\eta_{\rm
BBN}=(6.143\pm 0.190)\times 10^{-10}$ and $\eta_{\rm CMB}=(6.104\pm
0.058)\times 10^{-10}$ for 3 neutrino species). The latest \textit{Planck} best
fit value of the baryon density for the TT+TE+EE+lowE+BAO data combination is
$\omega_b=0.02242 \pm 0.00014$ \citep[Table\,2 of][]{Planck2018parameters}.

We define $\delta\eta/\eta=(\eta_{\rm BBN}-\eta_{\rm CMB})/\eta_{\rm CMB}$.
Then, after eq.~(\ref{eq:eta}) and the analogous equation at decoupling, the
fractional difference, $\Delta\epsilon/\epsilon_\gamma$, between the photon
energy density at decoupling ($\epsilon_{\gamma,\rm CMB}$) and at BBN
($\epsilon_{\gamma,\rm BBN}$), can be written as:
\begin{equation}
\frac{\Delta\epsilon}{\epsilon_\gamma} \equiv \frac{\epsilon_{\gamma,\rm CMB}-\epsilon_{\gamma,\rm BBN}}{\epsilon_{\gamma,\rm CMB}} \approx
1-\left(\frac{\eta_{\rm CMB}}{\eta_{\rm BBN}}\right)^{4/3}\simeq \frac{4}{3}\frac{\delta\eta}{\eta},
\end{equation}
where we have  assumed $(\Delta\epsilon/\epsilon_\gamma, \,\delta\eta/\eta) \ll
1$ and used the fact that $\epsilon_\gamma\propto n_\gamma^{4/3} \propto
\eta^{-4/3}$.

$\Delta\epsilon/\epsilon_\gamma$ is consistent with zero with an r.m.s.
uncertainty of
\begin{equation}
\sigma(\Delta\epsilon/\epsilon_\gamma)\simeq\frac{4}{3}\left[(\delta\ln\eta_{\rm BBN})^2+
(\delta\ln\eta_{\rm CMB})^2\right]^{1/2}\simeq 0.039,
\end{equation}
so that the 95\% confidence upper limit to the fractional amount of energy
density that could have been added to the CMB after the BBN is 0.078.

A different approach to derive model independent constraints on electromagnetic
energy releases after primordial nucleosynthesis was adopted by
\citet{Chluba2020}. Their argument goes as follows. Deep in the radiation
dominated era, but after the electron/positron annihilation, the expansion
timescale is controlled by the energy density of photons and neutrinos,
$\rho=\rho_\gamma + \rho_\nu$. A difference, $\delta\rho$, between the value of
$\rho$ at BBN and its standard value, due to a subsequent energy injection,
leads to a non-standard expansion rate, impacting on the primordial production
of light elements.

Such difference can be parameterized by a variation of the effective number of
neutrino species, $N_{\rm eff}$ \citep{SimhaSteigman2008}. For
$\delta\ln(N_{\rm eff})\ll 1$ \citet{Chluba2020} found $|\delta\rho/\rho|\mathrel{\rlap{\lower 3pt\hbox{$\sim$}}\raise 2.0pt\hbox{$<$}}
|0.2417\delta\ln(N_{\rm eff})/(f_\nu - 0.4089)|$, where $f_\nu$ is the fraction
of energy injected going into neutrinos. \textit{Planck} measurements of the
CMB anisotropy power spectra yielded $N_{\rm eff}=2.99\pm 0.17$ ($1\,\sigma$
error). Using the 95\% confidence uncertainty, $\delta\ln(N_{\rm eff})=0.112$,
\citet{Chluba2020} found $\Delta\epsilon/\epsilon_\gamma \mathrel{\rlap{\lower 3pt\hbox{$\sim$}}\raise 2.0pt\hbox{$<$}}
0.11$\footnote{The value $\Delta\epsilon/\epsilon_\gamma \mathrel{\rlap{\lower 3pt\hbox{$\sim$}}\raise 2.0pt\hbox{$<$}} 0.077$ quoted
in that paper is a misprint (Jens Chluba, private communication).} in the case
of a release of pure electromagnetic energy ($f_\nu=0$), and
$\Delta\epsilon/\epsilon_\gamma \mathrel{\rlap{\lower 3pt\hbox{$\sim$}}\raise 2.0pt\hbox{$<$}} 0.25$ for $f_\nu=0.5$, both limits being
at the 95\% confidence level. \citet{Chluba2020} also point out that if $f_\nu
\simeq 0.4089$  the energy release is unconstrained by measurements of $N_{\rm
eff}$.

A more precise calculation of constraints on $\Delta\epsilon/\epsilon_\gamma$
should take into account that they depend on the difference between the
effective number of neutrino species derived from BBN \citep[$2.86\pm
0.15$;][]{Fields2020} and that derived from CMB anisotropies, referring to the
much later decoupling epoch \citep{SimhaSteigman2008}. The error on the $N_{\rm
eff,BBN}-N_{\rm eff,dec}$ is a factor of 1.33 higher than the error on $N_{\rm
eff,dec}$ used by \citet{Chluba2020}. This also implies that the future
strengthening of the constraints on $\Delta\epsilon/\epsilon_\gamma$ depends on
decreasing the uncertainty on $N_{\rm eff}$ at both epochs.

The 95\% confidence limit on $\Delta\epsilon/\epsilon_\gamma$ derived in this
paper is illustrated by the black horizontal line in
Fig.~\ref{fig:constraints}. Also shown, for comparison, are the limits obtained
by \citet{Chluba2020} from the uncertainty on $N_{\rm eff}$, for two values of
$f_\nu$. The lines on the left part of the figure show the constraints implied
by upper limits on $\mu$-type distortions, again from \citet{Chluba2020}.

At $z\mathrel{\rlap{\lower 3pt\hbox{$\sim$}}\raise 2.0pt\hbox{$<$}} 4-5\times 10^6$ constraints from COBE-FIRAS supersede those from
BBN. At higher redshifts, thermalization of additional energy injected into the
CMB smooths out distortions, thus weakening or erasing constraints on the
thermal history of the universe. Even releases of very large amounts of energy,
$\Delta\epsilon/\epsilon_\gamma \sim 1$, occurring at $z\simeq 6.5-7.5\times
10^6$ wouldn't leave imprints detectable by COBE/FIRAS \citep{Fixsen2009} or by
a PIXIE-like experiment \citep{KogutFixsen2020}, orders of magnitude more
sensitive.

\section{Conclusions}\label{sect:conclusions}

The constraint from primordial nucleosynthesis, derived in this paper, implies
that no more than 7.8\% of the present-day CMB energy density could have been
released after primordial nucleosynthesis. Much stronger constraints are set by
COBE/FIRAS measurements of the  CMB spectrum, but these are limited to $z\mathrel{\rlap{\lower 3pt\hbox{$\sim$}}\raise 2.0pt\hbox{$<$}}
5\times10^6$ \citep{Chluba2020}, i.e. to redshifts about two orders of
magnitude lower than the nucleosynthesis redshift.

Comparable constraints on post-BBN energy releases have been derived by
\citet{Chluba2020} from measurements of the effective number of neutrino
species. As pointed out by these authors, such constraints are affected by
uncertainties on the fraction, $f_\nu$, of energy going to neutrinos. They are
weak for $f_\nu \approx 0.3-0.5$ and can be completely avoided for
$f_\nu\approx 0.409$, a value close to the typical amount of energy carried by
neutrinos in TeV scale dark matter annihilations.

Constraints on the contribution to $N_{\rm eff}$ by high energy neutrinos from
dark matter decay in the early universe have been obtained by
\citet{AcharyaKhatri2020}.

Tighter constraints on early energy releases from specific processes have been
inferred from measurements of the light element abundances. Examples are
Hawking evaporation of primordial black holes \citep[e.g.,][]{Keith2020} and
decaying particles \citep[e.g.,][]{Kawasaki2020}. However these constraints are
somewhat model dependent as they rest on the details of hadronic and
electromagnetic interactions involved. On the contrary, the limit presented
here is general and model independent.

Above $z\sim 5\times10^6$ thermalization effects wash out spectral distortions
induced by energy injections into the CMB. Planned next-generation experiments
like PIXIE \citep{Kogut2019, Chluba2019}, PRISM \citep{Andre2014},
PRISTINE\footnote{\url{https://www.ias.u-psud.fr/en/content/pristine}} or the
microwave spectro-polarimetry mission proposed by \citet{Delabrouille2019} will
reach sensitivities orders of magnitude higher. However the redshift range over
which spectral distortions are visible will be only marginally extended.

\section*{Acknowledgements}
We are grateful to the referee and to Jens Chluba for useful comments. Jens
Chluba also provided helpful clarifications on the \citet{Chluba2020} paper. MB
acknowledges support from INAF under PRIN SKA/CTA FORECaST and from the
Ministero degli Affari Esteri della Coo\-pe\-ra\-zione Internazionale -
Direzione Generale per la Promozione del Sistema Paese Progetto di Grande
Rilevanza ZA18GR02.

\section*{Data Availability}

The research described in this letter does not make use of any database. The
data used are those available in the references cited when the data are
introduced. The data generated are given in the text.

\bibliographystyle{mnras}