Visible to Mid-Infrared Optical Constants of Orthopyroxenes

Radiative transfer models of remotely acquired infrared spectra result in quantitative identification of minerals on planetary surfaces. Optical constants, or the real ( n ) and imaginary ( k ) indices of refraction are necessary inputs in such models. Pyroxenes are ubiquitous on the surfaces of terrestrial bodies within our solar system and can be readily used as thermo-barometers to interpret magmatic histories. However, optical constants for intermediate pyroxene compositions are undetermined. Here, we have determined the optical constants of two natural orthopyroxenes both in the visible/near-infrared (VNIR) and mid-infrared (MIR) regions. VNIR reflectance spectra were measured using powdered samples and modeled using a combination of Hapke theory and Kramers-Kronig analysis. MIR reflectance spectra were measured on oriented single crystal samples with non-normal incidence and modeled using Lorentz-Lorenz dispersion theory. The optical constants


Introduction
Pyroxenes are among the most abundant rock-forming minerals found in the Solar System and have been detected on the Moon, Mars, rocky asteroids, and interplanetary dust particles (IDPs).The structure of pyroxenes can be broadly subdivided into orthopyroxenes (OPX) and clinopyroxenes (CPX), having orthorhombic and monoclinic symmetry, respectively.Both structure types, consists of alternating tetrahedral and octahedral layers along the c-axis (Cameron & Papike, 1981).The octahedral layer includes two metal cation sites, M1 and M2.The M1 site allows for normal octahedral coordination, predominantly occupied by Mg 2+ and Fe 2+ , and the M2 site is in slightly disordered octahedral coordination, accommodating large cations, such as Ca 2+ (Morimoto et al., 1988).End members are: enstatite (En: MgSiO3), ferrosilite (Fs: FeSiO3), and wollastonite (Wo: CaSiO3).OPX are Mg-Fe rich, containing only up to 5% of wollastonite component.Monoclinic pyroxenes have higher contents of Ca.The distribution of iron between the M1 and M2 sites is a function of composition, temperature, pressure, and cooling history (Anovitz et al., 1988;Klima et al., 2007Klima et al., , 2008;;Molina & Zanazzi, 1991;Saxena et al., 1974;Virgo & Hafner, 1970).Thus, the quantitative characterization of pyroxene mineralogy of planetary surfaces is indicative of magmatic evolution and cooling conditions of the source rock (Klima et al., 2008;Lindsley, 1983).
The determination of pyroxene composition through electronic and vibrational spectroscopy, particularly in the near-to MIR (mid-infrared) range with applications to remote sensing, has been the subject of several studies (Adams, 1974;Cloutis & Gaffey, 1991a, 1991b;Hamilton, 2000;Klima et al., 2007Klima et al., , 2008;;Klima, Dyar, et al., 2011;Sunshine et al., 1990;Sunshine & Pieters, 1993).Observed IR spectral features result from both electronic and vibrational phenomena.The absorption, transmission, or reflection at a particular wavelength of light is specific to the ion and its bonding environment within a crystal lattice.Absorption in the visible to near-infrared (VNIR; ∼0.4-2.5 μm) is the result of electronic transitions involving Fe, whereas features in the MIR(∼2.5-50μm) are due to crystal lattice vibrations.
• Optical constants are necessary inputs into remote sensing models that quantitatively determine modal mineralogies on planetary surfaces • Many optical constants available are end member compositon only • Here we calculate the optical constants for two compositions of orthopyroxene, common minerals observed on planetary surfaces, in the visible/near infrared and mid to far infrared

Supporting Information:
Supporting Information may be found in the online version of this article.
Similar shifts are observed with increasing Ca 2+ content, where more drastic changes in the band 2 positions are seen for pyroxenes up to Ca 2+ saturation at the M2 site.Due to the primary occupancy of the M2 site by Ca 2+ , the Fe 2+ spectral trend can be masked for Ca-rich pyroxenes (>50% Wo; Klima, Dyar, et al., 2011).Several other cations affect the spectral signature of pyroxenes as well, including a shift of both major bands to shorter wavelengths with increasing aluminum content, for example, in terms of spectral shift, every 1 wt.% of Al 3+ is approximately equivalent to a 5% decrease in Fe 2+ (Cloutis & Gaffey, 1991b).Thus, VNIR spectroscopy alone is insufficient to quantitatively estimate pyroxene composition.
The MIR yields detailed structural information as the location and number of spectral features are directly related to the vibrational modes of repeating units within a crystal structure.In particular, the MIR is sensitive to Si-O vibrational modes (Arnold et al., 2014;Hamilton, 2000;Salisbury, 1972).This makes the MIR spectral range ideal for remote sensing applications Analysis in this region yields information about the whole composition of the material directly related to its structure.In the case of pyroxenes, spectral features are the result of Si-O vibrational modes as well as M-O vibrational modes (Farmer, 1974;Jäger et al., 1998).Analysis of MIR spectral data is complicated, however, by multiple scattering effects that result in reductions in fundamental band strength, the appearance of transparency features, and nonlinear spectral mixing when particle sizes are similar to or smaller than the wavelength of light (e.g., Salisbury & Wald, 1992).
Radiative transfer (RT) models are commonly used to analyze remotely sensed data to quantitatively determine the mineral compositions of planetary surfaces at wavelength ranges and grain sizes where linear mixing cannot be assumed.By directly modeling light interacting with the surface as a function of particle size, composition, roughness, incidence, and emission angles, etc., quantitative abundances can be estimated, provided that the necessary input parameters are known.Such models include those by Hapke (1981Hapke ( , 1996Hapke ( , 2012)), Lumme and Bowell (1981), Mishchenko et al. (1999), andShkuratov et al. (1999).An essential input into these light scattering models is the complex index of refraction, ñ = n + ik, where n and k are the real and imaginary indices of refraction, also known as optical constants.The optical constants are fundamental bulk properties of any given mineral that depend on both composition and crystal structure.
The two models most often used in planetary remote sensing applications are those of Hapke (1981Hapke ( , 1996Hapke ( , 2012) ) and Shkuratov et al. (1999).In Hapke's theory, n and k are used to calculate both the single scattering albedo (ω 0ν ) and the Ambartsumian-Chandrasekhar H-function, (H(x)), which are used to solve the radiative transfer equation to model reflectance (Chandrasekhar, 1960;Hapke, 1981Hapke, , 1996Hapke, , 2012)).In this case, while reflectance is not linear with composition under many conditions (specifically in the VNIR), single scattering albedo does vary linearly with percentage in an intimate mixture (either mol% or vol%, depending on several modeling factors), and thus the bulk single scattering albedo derived from a remote reflectance measurement can be used to extract quantitative abundance information.Skhuratov's model describes the reflectance (A) as a function of the four parameters n, k, S, and q, where n and k are the optical constants, S is the scattering path-length, and q is the particle packing density (Shkuratov et al., 1999).As noted by Sklute et al. (2015), k in Shkuratov's approximation of radiative transfer are dependent upon the path length light travels in the material.Therefore, the value calculated with this model, while potentially useful, is not strictly representative of the actual k of the material in question.Thus, while this model is useful in modeling surface composition if n and k are already known, it is not appropriate for the derivation of optical constants through inverse modeling.In either case, the successful use of RT models to make quantitative estimates of regolith surface composition and physical properties hinges on the availability of optical constants for a range of major rock forming minerals.
A summary of all pyroxene compositions for which optical data has been derived is shown in Table 1.In this work, the values of n and k have been determined using dispersion analysis for two unique OPX compositions in the range of 250-4,000 cm −1 (2.5-40 μm) following the methods outlined by Glotch et al. (2006Glotch et al. ( , 2007) ) and Arnold et al. (2014).In addition, the bulk VNIR optical constants for these same samples were calculated here following the procedure outlined in Sklute et al. (2015Sklute et al. ( , 2021) ) and Ye et al. (2021).Previous works deriving optical constants of OPX are summarized in Table 1.

Moon
The Clementine mission enabled multispectral VNIR global mapping of lunar surface composition, highlighting the heterogeneity of lunar regolith both laterally and with depth.Pyroxene was found to be the major component in mafic rocks, predominantly present in mare basalts and only as a minor component of the Mg-suite and ferroan anorthosite (FAN) rocks of the lunar highlands (Tompkins & Pieters, 1999).Cahill and Lucey (2007) and Cahill et al. (2009) used radiative transfer models of Clementine UVVIS and NIR data to study the global distribution of major minerals as a function of latitude, longitude, and crustal depth, finding exposures containing low-Ca pyroxene (LCP), at the South Pole-Aitken basin (SPA) and in a few crater peaks within the Procellarum KREEP (Potassium, Rare Earth Elements, and Phosphorus) Terrane (PKT).Expanding on these results (Klima, Pieters, et al., 2011), used the Modified Gaussian Model (MGM) to deconvolve NIR spectral data from the Chandrayaan-1 Moon Mineralogy Mapper (M 3 ) to constrain the distribution and composition of LCP with respect to candidate Mg-suite intrusions.LCPs were found to be regionally constrained to craters within SPA and to the north and south of Mare Frigoris in the northern hemisphere near side.Approximately half of the pyroxenes modeled are compositionally comparable to those found in conjunction with Fe-enriched (Mg # 55-75) FAN suite rocks.
The Diviner Lunar Radiometer Experiment aboard the Lunar Reconnaissance Orbiter (LRO) introduced complementary MIR multi-spectral data to the Clementine and M 3 VNIR spectral measurements, enabling the detection of iron-poor lithologies where mafic minerals are present.Diviner includes 3 channels devoted to mineral detection that have narrow passbands centered around 8 μm (Greenhagen et al., 2010).This region exhibits the emissivity maximum, or Christiansen feature (CF), associated with degree of SiO 2 polymerization where the CF positions shift to shorter wavelengths with an increase in SiO 2 polymerization (Conel, 1969;Greenhagen et al., 2010;Nash et al., 1993;Salisbury & Walter, 1989).Analysis of the CF region, as well as other major MIR features, including the Reststrahlen bands, is well suited to identify and distinguish Fe-poor lithologies from remote spectroscopic measurements.Diviner data have been used to develop global maps of silicate composition and space weathering across the lunar surface.While low spectral sampling has prevented definitive determination of mineral composition, CF positions of lunar mare material have been observed to be consistent with a pyroxene-dominated mineralogy (Greenhagen et al., 2010), with olivine (Arnold et al., 2016), nearly pure plagioclase (Donaldson Hanna et al., 2012), and other mineralogies (Glotch et al., 2010(Glotch et al., , 2011(Glotch et al., , 2015;;Lucey et al., 2021;Song et al., 2013) detected as well.
The upcoming Lunar Trailblazer mission, which aims to understand the Lunar water cycle, is to be launched with NASA's Interstellar Mapping Acceleration Probe (IMAP) in 2025, and includes two spectrometers, the High-resolution Volatiles and Minerals Moon Mapper (HVM 3 ) and the Lunar Thermal Mapper (LTM).The HVM 3 instrument has a spectral range of 0.6-3.6 μm and has been optimized to quantify water.The LTM includes four temperature channels and 11 compositional channels targeting the Si-O CF, which represents a significant increase in compositional channels compared to Diviner's three.Overall, these instruments will provide an infrared and MIR map datasets with an increase in both spatial and spectral resolution compared to previous missions (Ehlmann et al., 2021).

Mars
In addition to pyroxene identification in Martian meteorites (Papike et al., 2009), analysis of VNIR and MIR wavelengths showed LCP and HCP to dominate much of the mineralogy of the low-albedo regions of Mars.Global compositional mapping of Mars from the Mars Global Surveyor Thermal Emission Spectrometer (MGS-TES) identified two surface types in dark regions defined by modal mineralogy: Surface Type 1, mostly composed of basalt, containing ∼25% CPX based on linear mixing models; and Surface Type 2, basaltic andesite or weathered basalt, with ∼10% CPX content (Bandfield, 2000).Expanding on this work, Rogers and Christensen (2007) and Rogers et al. (2007) reexamined low-albedo regions using TES data and identified 4 surface types based on modal mineralogical differences from region to region, including refined definitions for Surface Type 1 and 2. The distributions found in these studies are consistent with VNIR observations made using the Phobos 2 Imaging Spectrometer for Mars (ISM) by the Mars Express Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité (OMEGA) VNIR data sets (Bibring, 2005;Bibring et al., 2006;Mustard et al., 1997).Both OPX and CPX are found in all 4 surface types (Rogers & Christensen, 2007).et al. (1970), interpreted the strong absorption band centered at 0.9 μm from telescope reflectance data (0.3-1.1 μm) of asteroid, four Vesta, as indicating a pyroxene-rich surface compositionally similar to basaltic achondrites.Subsequent studies found pyroxene to be a common surface component of S-and M-type asteroids (Chapman et al., 1975;Gaffey & McCord, 1978;Gaffey et al., 1993;Hardersen et al., 2005).Similarities between the spectral signature of Vesta and that of the howardite-eucrite-diogenite (HED) meteorites suggest Vesta to be the parent-body of HEDs (Drake, 1979;Feierberg & Drake, 1980;Gaffey, 1997).Recent support for this conclusion comes from spectral maps provided by the Visible and InfraRed (VIR) imaging spectrometer on the Dawn spacecraft.De Sanctis et al. ( 2012) provided geologic context for HED meteorites by highlighting the spectral differences with depth and across Vesta's surface with respect to HED compositions.The average spectral signature at global to regional scales is dominated by a howardite-like regolith containing localized regions with various proportions of eucrite (∼50% CPX) and diogenite (∼90% OPX) lithologies.VIR data of the Rhea Silvia Basin, a deeply excavated impact basin in the South Polar Region, indicates a larger amount of diogenite-like material at the surface suggesting diogenite may dominate at depth.Exogenous Vesta or Vesta-like material has also been detected on the near-Earth asteroid Bennu by the OSIRIS-REx OVIRS instrument, demonstrating the prevalence of asteroidal mixing processes in the inner solar system subsequent to planetary accretion (DellaGiustina et al., 2021).Sunshine et al. (2004) demonstrated that high ratios of HCP/(HCP + LCP) are indicative of a differentiated body, using VNIR data collected at the NASA Infrared Telescope Facility (IRTF) and analyzed with the MGM.This was applied to several asteroids within the Merxia and Agnia S-type asteroid families and Thetis (S-type, non-family), all of which were found to have high HCP/(HCP + LCP) ratios (0.27-0.64) indicating they were derived from differentiated bodies.

Interplanetary, Cometary, Circumstellar, and ISM Dust
Enstatite (Mg-rich OPX) was observed in interplanetary dust particles collected from the Earth's stratosphere (Bradley et al., 1983;Mackinnon & Rietmeijer, 1987;Sandford & Walker, 1985), and in comets through groundbased IR observations (Hanner et al., 1994).Apart from these exceptions, however, it was commonly maintained that cosmic dust silicates are amporphous.Thus, much of the analytical work focused on deriving optical constants for glasses of pyroxene composition (Dorschner et al., 1995;Henning & Mutschke, 1997;Jäger et al., 1994).However, when the high-resolution MIR and FIR spectra from the Infrared Space Observatory (ISO) in the 2.5-240 μm wavelength range became available, crystalline pyroxenes were identified around young and evolved stars (Malfait et al., 1998;Molster et al., 2002), in the circumstellar disks of young stellar objects (Bouwman et al., 2001), and in planetary nebulae (Waters et al., 1998).With these new observations, it became apparent that there were insufficient high-resolution laboratory spectra of Mg-Fe silicates, against which to quantitatively interpret ISO data.Several studies were therefore undertaken to remedy this and understand the spectral properties of silicates, including olivine and pyroxene, across the ISO spectral range (Chihara et al., 2002;Jäger et al., 1998;Koike et al., 2000).

Sample Description and Preparation
Two large natural single crystal samples, enstatite (En1; Bamble, Norway) and ferroan-enstatite (Fen1; Labrador, Canada), were used for this study with an approximate size of 0.35-1 cm on a side (Table 2).Samples were loaned to us from the mineral collection at Stony Brook University, Department of Geosciences.Chemical compositions were determined by the Cameca SX100 electron microprobe at the American Museum of Natural History, New York.Thin sections of each sample were prepared to optically identify any possible exsolution or twinning that might affect IR data.While neither of these two samples exhibit any apparent twinning or exsolution, sample En1 does have parting parallel to the 001 plane.The IR spectra for En1 agree with oriented reflectance spectra for a natural enstatite (En 96 Fs 4 ) reported by Jäger et al. (1998) and synthetic ortho-enstatite reported by Demichelis et al. (2012), indicating that the parting did not affect our spectral results.
Each sample was cut using a diamond saw blade to produce three different surfaces parallel to the a-plane (100), b-plane (010), and c-plane (001).Each surface was polished to a 0.25 μm surface roughness to provide a mirror-like surface.The orientations of these surfaces were confirmed through single-crystal X-ray diffraction measurements from which each crystal face was indexed.This analysis consisted of determining the unit cell and orientation matrix by the collection of a preliminary set of frames.From this, we indexed each crystal face using the CrysAlisPRO software.The reflections were collected using a four-circle κ Oxford Gemini diffractometer with an Atlas detector (λ = 0.71073Å) in the Department of Chemistry at Stony Brook University.
The sieving process was repeated three or more times.Each grain size was washed in DI water to remove clinging fines from larger grains.Powder X-ray diffraction patterns were collected to assess sample purity using a Miniflex diffractometer with CuKα radiation (λ = 1.54059Å) and angle dispersive geometry, where data was collected in the 2θ 3-100°, with a step size of 0.02°, and continuous collection.A pressed pellet of 13 mm diameter and ∼1.5 mm thickness was made from the smallest size fraction of each sample.

Mid-and Far-IR Spectra
We collected MIR specular reflectance spectra of oriented single crystals using the Stony Brook University Center for Planetary Exploration's Nicolet 6700 Fourier transform infrared (FTIR) spectrometer equipped with a FT-30 specular reflectance accessory with 30° incidence and emergence angles at a spectral sampling of 2 cm −1 .MIR spectra were measured using a KBr beamsplitter and a DTGS detector with a CsI window in the 400 to 4,000 cm −1 wavenumber range.FIR spectra were collected using a Nicolet Solid Substrate beam splitter and a DTGS detector with a polyethylene window in the 250 to 600 cm −1 wavenumber range.A total of 256 and 512 scans were averaged to produce the MIR and FIR spectra respectively.A KRS-5 holographic wire grid polarizer was placed in the path of the light source beam for the collection of polarized IR reflectance spectra for the single crystal samples.During analysis, the samples were placed on an aperture mask that controlled the spot size.The mask spot size was selected as the largest to be completely covered by the polished sample face.The two aperture sizes are approximately 4.7 and 6.8 mm and were used for Fen1 and En1 respectively.The underside of each mask was coated with paraffin soot prior to data collection to reduce light scattering from the mask.A background spectrum was acquired using a first surface gold mirror.A "blank" spectrum was also acquired of the empty mask and was subtracted from the mineral spectrum.The signal produced by the mask did not have a major effect in the MIR range, but a steep incline was observed beginning around 470 cm −1 increasing the reflectance from ∼0.01 to ∼0.11, as shown in Figure S1 in Supporting Information S1.We collected a total of 3 spectra per sample with the polarization direction parallel to each crystal axis (E||a, E||b, and E||c).The MIR and FIR spectra for each orientation were joined between 600 and 700 cm −1 , and the overlapping region was used to scale the FIR spectral contrast to be consistent with that of the MIR spectra.
The smallest grain size fraction for each sample was pressed into a powder pellet ∼1.5 mm thick and 13 mm in diameter.Unpolarized MIR reflectance spectra of these pellets were collected using the same setup described for the MIR spectra of the oriented single crystal samples.

VNIR Spectra
VNIR reflectance spectra of the three largest size-fraction powders were measured using an ASD FieldSpec3 Max spectroradiometer using a 512 element Si photodiode array detector for the 0.35-1 μm (∼28,571-10,000 cm −1 ) range and two TE cooled InGaAs photodiode detectors for the 1-2.5 μm (∼10,000-4,000 cm −1 ) range.The samples were loaded into an aluminum sample cup painted flat black.The bottom of the sample cup was gently tapped on a table to level the powder's surface.Spectra of each size fraction for both samples were collected at 7 phase angles every 15° from 15° to 105° in the absence of ambient light, illuminated by a 50 W quartz halogen lamp, referenced to a calibrated Spectralon standard, and corrected for the non-Lambertian and wavelength-dependent nature of Spectralon following the methods of Yang et al. (2019).To reduce systematic error, the incident light source was kept at 45° to the sample, and the fiber optic bundle collecting the reflected light was adjusted on a goniometer to produce the desired phase angle.A total of four spectra were collected at each angle by turning the cup 90° for each collection to reduce error associated with possible preferred orientation of grains or uneven surface.The average of these four spectra was used in the model analysis discussed in Section 3.4.

Modeling of Mid-IR Oriented Optical Constants
MIR optical constants were derived using a combination of Lorentz-Lorenz dispersion theory and a Fresnel reflectance model (Glotch et al., 2007).Dispersion theory mathematically defines the vibrational transitions within a mineral as the summation of lattice harmonic oscillators.Each oscillator is described by a set of parameters, ν, 4πρ, and γ, which are the center frequency of oscillation, the band strength, and the band width, respectively.An additional term, the high frequency dielectric constant (ε 0 ), is a bulk property and can be approximated as   2 vis (Roush et al., 1991).The real (n) and imaginary (k) indices of refraction can be described as a function of these parameters by the following equations: The optical constants are then related to the total reflectance by the Fresnel equations, which describe the total reflectance for non-normal incidence angles (Glotch et al., 2007): (3) We used Matlab non-linear least squares optimization routine to iteratively fit measured spectra.Initial guesses for the parameters ε 0 , ν, 4πρ, and γ, were estimated through visual inspection of spectral features and input into the model.These values were then used in the model to iteratively fit the measured spectrum and calculate n and k values, which were used to calculate a modeled reflectance spectrum using the Fresnel relationships shown above in Equations 3 through 5.The model output includes optimized values of the input parameters, a calculated reflectance spectrum, n, and k.The optical constants are derived when a reasonable fit (<10% difference from model) of the measured spectrum is achieved using the least number of oscillators possible (Spitzer & Kleinman, 1961).

Modeling of VNIR Optical Constants
VNIR optical constants were derived following a modified procedure from Sklute et al. (2015Sklute et al. ( , 2021)).This model utilizes the Hapke treatment of radiative transfer combined with Kramers-Kronig transformations to iteratively fit measured spectra and obtain bulk values for n and k over the VNIR wavelength range (0.35-2.5 μm).
The reflectance is calculated from: where REFF is the reflectance of the sample ratioed to the reflectance of a Lambertian scatterer viewed at the same angle, K is the porosity factor, μ is the cosine of the emission angle, μ 0 is the cosine of the incidence angle, w is the single scattering albedo, p(g) is the phase function, and H(x) is the Ambartsumian-Chandrasekhar's H-function.
The phase function, p(g), represents the phase angle dependence of singly scattered light and is modeled with a two term Legendre polynomial: where g is the phase angle (the sum of the incidence and emission angles), and b and c are the phase function coefficients.
The measured radiance is related to n and k through the single scattering albedo, w, which is described as follows for closely packed particles: where Θ, S e and S i are the internal transmission factor and Fresnel reflectance coefficients, respectively.The internal transmission factor, Θ, takes into account the effective grain size D, that is, the path length traveled by light through a grain, the absorption coefficient (α = 4πk/λ), and the internal scattering factor, s, such that, The Fresnel reflectance coefficients (Equation 10) are integrated over all angles and can be described as a function of n and k as follows: The Ambartsumian-Chandrasekhar's H-function can be approximated by or its exact solution can be determined using a lookup table approach.For highly reflective materials, the exact solution is necessary for good convergence.For the samples in this study, the approximation (Equation 15) was used.
The greatest challenge to modeling the VNIR optical constants is that without additional constraints, one can derive a model (Equation 8) which appears to approximate the data well using multiple mathematically equivalent solutions for variables n(λ), k(λ), b, c, D, and s, providing non-unique solutions, many of which might not be physically reasonable.Following Sklute et al. (2015), we used spectra of three grain sizes at 7 phase angles each to derive the intrinsically wavelength and grain size independent optical constants n and k.However, we reordered the steps in a manner similar to Yang et al. (2019) as outlined in Sklute et al. (2021) to ensure that phase behavior is completely described by phase function coefficients, thus obtaining a more robust solution for k.The procedure is as follows: 1.Because each grain size has one single scattering albedo, w, at each wavelength, all phase angle spectra for each grain size can be modeled using the same w, thus forcing all phase behavior to be accounted for by the phase function coefficients b and c.Once wavelength dependent b and c are determined for each grains size through this step, they are considered known quantities for the remainder of the derivation.2. The single scattering albedos from the three grain sizes viewed at a single-phase angle (typically g = 30°, although solutions from all phase angles at this stage are compared to look for systematic variations), along with a constant n obtained from the average value n at the sodium D line, are used to derive a wavelength-dependent value of k.We thus force all variation between the spectra of different grain sizes to be accounted for by effective grain size and internal scattering coefficient (Equations 10-15).3. The derived VNIR k is combined with the MIR k and a Singly Subtractive Kramers-Kronig (SSKK) transformation is used to derive the wavelength-dependent VNIR n: While the SSKK transformation ideally needs an infinite data range to converge to a unique solution, this condition is never met.Sklute et al. (2015) demonstrated that extending the data range in the MIR for calculation of VNIR optical constants results in substantially different n values than if only the VNIR range is used.
Here, we extend the k data to the MIR by modeling the pressed pellet spectrum (effectively averaging over all possible orientations) in the same manner as the oriented spectra.
4. The wavelength dependent n is used in the minimization procedure to derive a new k and these steps are iterated several times until there is negligible variation between iterations.

Results
Here we provide the observed and modeled near-IR and mid-to far-IR reflectance, as well as the calculated values of n and k for both samples.All spectra and derived optical constants are available at https://doi.org/10.5281/zenodo.4758045

Mid-IR Oriented Optical Constants
The oriented observed and modeled reflectance spectra for enstatite (En1) are shown in Figure 1  enstatite at all three orientations are shown in Table 4.The fit for the ferroenstatite MIR E || b has some misfit at low wavenumbers.The spectrum was best fit with an oscillator at the end of the available data: there may be an oscillator here, but without an extended spectrum this in unconfirmed and is reported here with caution.

VNIR Optical Constants
Observed and modeled reflectance for all three grain sizes of En1 in the VNIR are shown in Figure 5.For clarity, only fits for a phase angle of 45 are shown.Fits for all other phase angles are shown in Figure S2 in Supporting Information S1.The derived optical constants, n and k, are shown in Figure 6 Observed and modeled reflectance for all three grain sizes of Fen1 in the VNIR are shown in Figure 7, for clarity, only fits for a phase angle of 45 are shown.Fits for all other phase angles are shown in Figure S3 in Supporting Information S1.The derived optical constants, n and k, are shown in Figure 8.

MIR
The typical approach taken for deriving optical constants in the MIR, particularly in the case of poorly crystalline materials, that is, clays and iron oxides (Glotch & Rossman, 2009;Roush et al., 1991), involves deriving n and k from reflectance spectra of pressed pellets.These optical constants were assumed to be the average of all principal indices of refraction.While this is a good approximation, it is not completely accurate, as can be seen in Figure 9 where we show the average of all three oriented indices of refraction calculated here for En1 set against the calculated n and k from a pressed pellet of the same sample.Specifically, n and k values for the pressed pellet tend to be lower than the average of the three orientations of the single crystal by a factor of ∼2 in the Reststrahlen band region.From the standpoint of deriving VNIR optical constants, the most important feature of the MIR k is its magnitude near 2.5 μm, as this has the greatest effect on the degree to which VNIR n decreases with λ.

Note.
Oscillators are listed for EIIa, EIIb, EIIC separately.The standard error (σ) is given with each parameter.

Table 3
Oscillator Parameters for Enstatite

Note.
Oscillators are listed for EIIa, EIIb, EIIC separately.The standard error (σ) is given with each parameter.

Table 4
Oscillator Parameters for Ferroenstatite

10.1029/2021EA002104
14 of 23 There has only been two previously reported oriented optical constants derived for orthopyroxene previously reported in the literature (Jäger et al., 1998;Zeidler et al., 2015).( Zeidler et al., 2015) investigated the temperature dependance of oriented optical constants of orthoenstatite (En 92 ) using a Lorentz oscillator model over a spectral range of 5-60 μm (∼150-2,000 cm −1 ).The oscillator ν j positions presented here are in good agreement with those reported in Zeidler et al. (2015), although some shifting of this parameter is present.This may be due to compositional effects that are expected to shift the oscillator positions systematically with increasing iron content as observed in bulk emissivity data (Hamilton, 2000).While beyond the scope of this study, the compositional effect to individual oscillator parameters is of interest.

VNIR
Other studies have derived k values for OPX by converting NIR diffuse reflectance data to k spectra using different relations (Lucey, 1998) or converting spectra to k using MGM (Denevi et al., 2007;Trang et al., 2013).There are four key differences between the methodology used here and those that have been previously reported.First, in other studies, the k values are derived for each grain size individually and averaged together.By simultaneously fitting our multiple grain size data to derive a single k spectrum, we arrive at a more accurate solution for a grain size independent imaginary index of refraction.Second, the wavelength-dependent n is also calculated here, whereas previous studies have assumed this to be constant.While the change in n across this spectral region is small, we have shown here there is some variation of n with wavelength (Figure 6).Third, in our calculations of VNIR optical constants, we extended the data range to the MIR to reach a more convergent and correct solution for the SSKK transformation (Sklute et al., 2015).Finally, and most importantly, we have leveraged the fact that grain size contributions and phase behavior contributions are mathematically separable using Hapke theory in order to separate phase function determination from optical constant determination.Other studies assumed either that phase behavior is isotropic or have chosen constant values for the phase function coefficients, which forces any variation due to phase behavior to be accounted for by the optical constants themselves.Thus, our calculated n and k in the VNIR are more robust than previous estimates.

Mid-IR
Following Arnold et al. (2014), we assessed the MIR dispersion model fits using standard error for each of the oscillator parameters and reported the errors in Tables 3 and 4 for En1 and Fen1 respectively.The inherent error associated with the parameters used to calculate the final n and k values was estimated by shifting the parameter values by ±20 cm −1 for ν j and ±50% for others.Little to no change is seen in the resulting n and k values.Representative spectra of the resulting n and k for En1 E||a with the parameter changes are shown in Figures S4-S7 in Supporting Information S1.With each change there are only minor differences that are associated with smaller or clumped spectral features.

VNIR
The assumptions we made in our approach to simplify the manner in which we mathematically describe these material properties allowed us to achieve a great fit for the VNIR reflectance data, shown in Figures 5 and 7 and Figures S2 and S3 in Supporting Information S1, adding confidence in the calculated k values.Minimal misfit is observed in the UV region, see Figures S2 and S3 in Supporting Information S1, likely due to the extension of data into this region through a linear fit rather than direct measurement.
In order to assess the accuracy of the model, k values were calculated for each phase angle individually (see step 2 in model outline) and compared in Figure 11.There are two major features in the spectra at ∼0.9 and ∼1.9 microns.The positions of the maxima do not shift, however the k values vary with phase angle, across the entire spectral range, with the greatest shift observed at the maxima.According to Equations 10-14, k values are independent of phase angle, thus no variation should be observed.The difference in intensity in Figure 11 seems to be loosely correlated with phase angle, where the smallest and largest phase angles have lower calculated k values than for those that are in the mid-range of phase angles used here.We observe inverse behavior in the measured reflectance spectra.Figure 12 shows a comparison of k value and the reflectance for all grain sizes at the ∼0.9 μm maximum for both En1 and Fen1, where a pseudo-parabolic relationship between the phase angle and reflectance can be seen (additional data are shown in Figure S8 in Supporting Information S1).
A possible explanation for this behavior is that the Spectralon is an isotropic scatterer, rather than a Lambertian scatterer (Sklute et al., 2015).An isotropic scatterer would have a higher reflectance at mid phase angles and a lower reflectance at high and low phase angles with respect to a Lambertian scatterer (Piatek, 2003;Sklute et al., 2015).The model would, therefore, attribute more reflectance to samples with respect to the standard at mid phase angles.This would result in a systematic difference in the value of k at all wavelengths, which is what we observe here for both En1 and Fen1, however may not account completely for the greater degree of spread in k value at the feature maxima.
Additionally, both the single scattering albedo and the k solution are influenced by the calculated grain size values, which may contribute to this affect.These values are calculated within our model and are allowed to vary within lower and upper bounds based on the range of grain sizes for each sample.The resulting calculated grain sizes for each phase angle are shown in Table 5.In the case of enstatite, the best fit is achieved when both the small and medium grain size samples calculated to have similar mean grain sizes of ∼63 μm.The calculated grain size for the smalles grain size sample pins at the top of the grain size range (63 μm), regardless of phase angle.This may be due to over-processing of the sample, where too much of the smaller grains were removed during the fines removal process, skewing the overall grain size to higher than expected.However, there are other explanations.Grain size is one of the largest contributing factors to overall reflectance (Helfenstein & Shepard, 2011).If some other phenomena were darkening the spectrum, the model would likely use grain size to reflect that change.This could include small dark inclusions, admixing of the sample cup from sample of insufficient depth or spot size, among others.In the end, it is important to remember that this parameter is related to grain size but is mathematically the average path length traveled by the light before being scattered, reflected, or absorbed.The large grain size of enstatite fits well within the expected range with a minimal variation with phase angle (Table 5).In the case of ferroenstatite, there is a clear relationship between larger overall k values and lower grain size estimates, however there is not a deviation of overall reflectance with phase angle that could account for this behavior.Ultimately, the model should arrive at a single grain size value regardless of phase angle, but this relationship between larger grain size estimates and lower reflectance values is likely influenced by the isotropic nature of the Spectralon as well as other phenomena not explicitly modeled here, such as microscopic surface roughness or shadow hiding.An in-depth study of the effect of grain size, grain shape, sample packing, and sample composition is beyond the scope of this paper but would shed light on these modeling results.Another contribution to model errors we must consider in our approach is the porosity correction.As discussed by Sklute et al. (2015), there is a large degree of uncertainty with its value if accounted for in the model.For this reason, we report n and k with an uncorrected porosity factor.However, this can result in the k values being too   small by up to a factor of 2 (Hapke, 2002).Using the raw data provided with this paper along with the computer program available as part of this work, the interested reader is able to explore the interaction between n, k, and porosity.The optical constants presented here were calculated using a set of assumptions common in the literature (e.g., Roush et al., 1991).Decreasing the number of assumptions necessarily increases the number of modeling parameters, thus affecting the uniqueness of the solution.The methodology used in this paper seeks to strike a balance between limiting assumptions and uniqueness of fit while providing the tools and data for further exploration.

Conclusions
We have provided a comprehensive data set of optical constants of two compositions of natural OPX for use by the Earth and planetary remote sensing communities.This study focused on two distinct compositions and, while other studies have concentrated on the NIR with synthetic ortho-and clino-pyroxenes across a wide range of compositions (Denevi et al., 2007;Lucey, 1998;Trang et al., 2013), there is still a lack of oriented optical constant data in the MIR for more iron-rich orthopyroxene species.This is mostly due to the lack of natural or synthetic single crystals of Fe-rich OPX large enough for current benchtop analytical capabilities.With the detection of Fe-rich OPX in eucritic meteorites (Pang et al., 2016), it is possible this mineral species may be an important component on other differentiated asteroids.Thus, future work widening the available orthopyroxene data set to include oriented MIR data for Fe-rich compositions will aid in these studies.Additionally, the anomalous behavior observed here in the variation of VNIR k with phase angle at the absorption maxima could not be completely explained as no systematic effects of the modeled reflectance or grain size estimates could be directly correlated with the observed behavior.While tracking the source of this discrepancy is outside of the scope of the current work, future work should be conducted to explain these aspects of Hapke theory.The updated Sklute et al., 2021 Matlab code is available for download at https://zenodo.org/record/4429127.
, where the incident polarization directions are parallel to the axes of symmetry of the crystal, a-(top), b-(middle), and c-(bottom), denoted by E || a, E || b, and E || c respectively.The derived n and k values for all three polarization directions are shown in Figure 2. Oscillator parameters for enstatite at all three orientations are shown in Table 3.The oriented observed and modeled reflectance spectra for ferroenstatite (Fen1) are shown in Figure 3.The derived n and k values for all three polarization directions are shown in Figure 4. Oscillator parameters for
Jäger et al. (1998) used a Kramers-Kronig analysis to calculate the n and k optical constant values for a natural single crystal of enstatite (En 96 Fe 4 Wo 0 ).In Figure10, we show a comparison of the n and k for enstatite in this study (En1) and that reported in Jäger

Figure 6 .
Figure 6.Average Visible to Near-Infrared optical constants, n and k, of enstatite (En1).Separate n and k values were calculated for every phase angle and averaged to find the reported value.

Figure 8 .
Figure 8.Average Visible to Near-Infrared optical constants, n and k, of ferroenstatite (Fen1).Separate n and k values were calculated for every phase angle and averaged to find the reported value.

Figure 9 .
Figure 9.Comparison between enstatite (En1) optical constants, n and k, from the average of the oriented Mid-IR values (blue, solid), and values calculated using a pressed powder pellet of the same sample (green, dotted).

Figure 10 .
Figure 10.Comparison of enstatite optical constant, n and k, with E-light polarization parallel to the c-axis (blue, solid) with optical constants reported for oriented enstatite of similar composition (green dotted), reported in Jäger et al. (1998).

Figure 11 .
Figure 11.Individual calculated k values at each phase angle for enstatite (top) and ferroenstatie (bottom).

Figure 12 .
Figure12.Maximum k value (filled circles) at the ∼0.9 μm with phase angle for enstatite (top) and ferroenstatite (bottom).Maximum reflectance value for the small grain size (triangles), medium grain sizes (diamonds) and large grain size (pentagon) at the ∼0.9 μm for enstatite (top) and ferroenstatite (bottom).Filled reflectance symbols represent measured reflectance, and hollow black symbols represent the modeled value.