Simulation of High-contrast Polarimetric Observations of Debris Disks with the Roman Coronagraph Instrument

The Nancy Grace Roman Space Telescope Coronagraph Instrument will enable the polarimetric imaging of debris disks and inner dust belts in the optical and near-infrared wavelengths, in addition to the high-contrast polarimetric imaging and spectroscopy of exoplanets. The Coronagraph uses two Wollaston prisms to produce four orthogonally polarized images and is expected to measure the polarization fraction with measurement errors <3% per spatial resolution element. To simulate the polarization observations through the Hybrid Lyot Coronagraph (HLC) and Shaped Pupil Coronagraph (SPC), we model disk scattering, the coronagraphic point-response function, detector noise, speckles, jitter, and instrumental polarization and calculate the Stokes parameters. To illustrate the potential for discovery and a better understanding of known systems with both the HLC and SPC modes, we model the debris disks around Epsilon Eridani and HR 4796A, respectively. For Epsilon Eridani, using astrosilicates with 0.37 ± 0.01 as the peak input polarization fraction in one resolution element, we recover the peak disk polarization fraction of 0.33 ± 0.01. Similarly, for HR 4796A, for a peak input polarization fraction of 0.92 ± 0.01, we obtain the peak output polarization fraction as 0.80 ± 0.03. The Coronagraph design meets the required precision, and forward modeling is needed to accurately estimate the polarization fraction.


INTRODUCTION
Despite recent advances in high-contrast imaging, circumstellar debris disks around main sequence stars are still poorly understood.Debris disks are composed of planetesimals, predominantly of dust with a small percentage of gas, resulting from successfully formed planetary systems.The analysis and study of debris disks provide valuable insights into the planet formation process and the structure of a planetary system (Backman 2004;Hughes et al. 2018).Our solar system can be broken into the inner hot and warm zodiacal dust, the cool asteroid belt, and the Kuiper belt that form the debris disk of our solar system (Wyatt & Jackson 2018;Levasseur-Regourd et al. 2020).
The boundaries of these populations are sculpted by the gravitation influence of solar system planets, which suggests an indirect technique for detecting planets via gap clearing (Stark & Kuchner 2008;Kennedy & Piette 2015), a distinct process from planet-driven gas shocks which open gaps in protoplanetary disks (Bae et al. 2017).The size of observed gaps in debris disks will depend heavily on dust composition since, for a fixed planet mass, the size of the gap depends on the transport rates, which is a function of the ratio of radiation pressure to gravitational attraction, commonly known as β.Accurately calculating β requires detailed knowledge of the dust grain properties: size, shape, composition, and porosity.In addition to the gaps, planetary companions are known to induce other features in the debris disks, such as warps, clumps, spirals, and brightness asymmetry (Wyatt 2008).Resolved, multi-wavelength observations of debris disks reveal general complementary information about the composition and morphology: optical and near-infrared (NIR) observations reveal scattered light by sub-micron and micron-sized dust grains, while observations in infrared (IR) and radio show thermal emission by sub-mm dust grains.
Although debris disks have been observed around a hundred stars in the last two decades, constraining disk properties has not been straightforward.Dust grains' properties contribute in a complicated and often degenerate way to radiative transfer processes, making it quite challenging to disentangle and constrain them individually using radiative transfer modeling (Krivov 2010).Scattered light from debris disks is expected to be linearly polarized due to asymmetries in their structure or scattering/absorption by the dust grains, where the polarization fraction (p = Q 2 + U 2 /I, where I, Q, and U are the Stokes parameters, (Stokes 1852)) as a function of scattering angle (the angle between the incident wave and the direction of the scattered wave) is also sensitive to specific dust grain properties such as composition, size, and distribution.Thus, when combined, total and polarized intensity measurements help constrain the geometrical and scattering properties of the debris disk more than with either one alone (e.g.Arriaga et al. 2020).Additionally, polarimetry can improve sensitivity to polarized sources relative to unpolarized star-light, improving the effective contrast ratio (Perrin et al. 2015).
Polarimetric observations of debris disks have been carried out using the Advanced Camera for Surveys coronagraph (ACS) in Hubble Space Telescope (HST) and current ground-based high-contrast imaging polarimeters in optical and NIR wavelength regions (Graham et al. 2007;Maness et al. 2009;Engler et al. 2017;Milli et al. 2017a;Esposito et al. 2020;Chen et al. 2020a;Hom et al. 2020;Crotts et al. 2021;Hull et al. 2022).Using the Gemini Planet Imager (GPI) (Macintosh et al. 2014) at the Gemini-South Telescope, Esposito et al. (2020) conducted a four-year survey of 104 stars and obtained polarization observations of 35 debris disks at NIR wavelengths.In addition to this, the Spectro Polarimetric High-Contrast Exoplanet REsearch (SPHERE)/Zurich Imaging Polarimeter (ZIMPOL) (Beuzit 2013;Schmid et al. 2018a), the SPHERE/The infrared dual-band imager and spectrograph (IRDIS) (Vigan et al. 2010), and the Nasmyth Adaptive Optics System -Near Infrared Imager and Spectrograph (NaCo) at the VLT (Witzel et al. 2011), and the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) (Martinache et al. 2016)/High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) (Hodapp et al. 2008) at the Subaru telescope have also been used to image debris disks in polarization at NIR and optical wavelengths to constrain disk properties (Milli et al. 2017a;Engler et al. 2017;Asensio-Torres et al. 2016).For example, polarimetric observations of HR 4796A using Gemini/GPI and VLT/SPHERE/ZIMPOL have helped to place some constraints on its dust grain composition (Arriaga et al. 2020;Milli et al. 2019).Similarly, using the total and polarized intensity observations through HST/ACS, Graham et al. (2007) identified distinctions between two approaches of modeling scattering phase functions, Henyey-Greenstein and Mie theory, to model the properties of the debris disk AU Mic.Although polarized observations have proven extremely useful in deriving the disk properties, it is quite challenging to simultaneously obtain the polarized and (unbiased) total intensity observations (Esposito et al. 2018).
The upcoming Nancy Grace Roman Space Telescope Coronagraph Instrument (Poberezhskiy et al. 2022(Poberezhskiy et al. , 2021;;Kasdin et al. 2020) will facilitate polarimetric observations of debris disks (sensitive to star-planet flux ratios of ≲ 10 −8 ) around nearby stars in addition to the high-contrast and highresolution imaging of exoplanets.The polarimetric module consists of two Wollaston prisms, each producing two orthogonally polarized images (I 0 , I 90 and I 45 , I 135 ) that are separated by 7.5".Polarimetric imaging is available for the narrow field with the Hybrid Lyot Coronagraph (HLC) at 575 nm and the wide field with the Shared Pupil Coronagraph (SPC) at 825 nm.The accuracy requirement in a linear polarization fraction (LPF) measurement per spatial resolution element (2×2 pixels in HLC and 3×3 pixels in SPC) is < 3%.However, Monte Carlo simulations using uncertainty in calibration, flat fielding, and photometric noise on standards estimate an RMS error in LPF measurement per resolution element to be 1.66% (Mennesson et al. 2021;Zellem et al. 2022).
Polarization observations through the Roman Space Telescope will play a vital role in providing constraints on some of the disks already observed by GPI and SPHERE with its complementary observations in optical wavelength regimes in addition to resolving the fainter dust rings much closer (around 1 AU) to nearby stars.One of the potential problems in polarimetry is the errors due to instrumental polarization and crosstalk arising due to the telescope and instrumental optics, which has been very well observed in GPI (Millar-Blanchaer et al. 2016;Millar-Blanchaer et al. 2022), SPHERE/IRDIS (de Boer et al. 2020;van Holstein et al. 2020), and SPHERE/ZIMPOL (Schmid et al. 2018b).To best plan the calibration strategies for Roman-CGI polarimetry observations, it is crucial to characterize the effect of instrumental polarization and crosstalk.Additionally, accurate end-to-end disk polarization observation simulations enable optimized observation planning through Roman-CGI.
In this work, we describe a process for simulating the polarization observations of debris disks through the Roman Coronagraph Instrument and demonstrate the instrument's potential to measure linear polarization fraction greater or equal to 0.3 with an uncertainty of less than 0.03.To demonstrate our approach, we generate simulated polarimetric observations of debris disk models around a nearby (early) Sun-like star Epsilon-Eridani (ϵ Eridani) and HR 4796A, an A0V star harboring an extensively studied debris disk 1 .As the main motivation of this paper is to demonstrate the simulations of disks through Roman-CGI, we use disk models that are analogues to the disks around these stars but do not necessarily agree with all the existing multi-wavelength observations.The paper is organized as follows: §2 describes the mathematical model for the observation simulation.The radiative transfer modeling approach of the ϵ Eridani and HR 4796A disks is described in §3.Generation of Point Response functions for the HPC and SPC mode is described in §4.§5 shows the creation of raw EMCCD images.The processing of disks incorporating noise and uncertainty factors from the observing scenario (OS) simulations is described in §6, and estimations of output polarization fractions are provided in §7.Finally, we provide our discussions and conclusions in §8.

MODELING METHODS
To give the reader tangible examples, we will describe the modeling process for two example systems, one so far unresolved in scattered light and one well-known, extensively studied.Figure 1 describes the outline of the different steps involved in the simulation of a polarization observation of a debris disk.First, we model the debris disk using the radiative transfer modeling software MCFOST (Pinte et al. 2006(Pinte et al. , 2009) ) to obtain the total intensity image and Stokes parameters Q and U images for linear polarization.Next, the orthogonal polarization components are convolved with the Roman Coronagraph point response functions (PRFs) obtained using the PROPER (Krist 2007) models run for the HLC and SPC modes.We use PROPER, which combines coronagraphic modes with EMCCD properties directly, instead of higher-level coronagraph models such as FALCO (Riggs et al. 2018) or CGISim.FALCO is mainly used for wavefront sensing and control, and CGISim is not required as we incorporate the EMCCD noise after the convolution.Next, the EMCCD raw images, including the EMCCD gain and noise characteristics, are generated.We then add the speckle and jitter noise using Observing Scenario 9 ("OS9")2 simulations of the Roman Coronagraph for the HLC mode and OS 11 ("OS11")3 simulations for the SPC mode.These images are processed as conventional CCD images with high read noise (no photon counting) using the "analog" mode as described in Nemati (2020).The final step is the estimation of Stokes parameters, polarization intensity, and total intensity after incorporating the instrumental polarization and polarization crosstalk from the pupil-averaged Mueller matrices4 of the instrument.Each step in the simulation is explained in detail in the following subsections.
Figure 1.The process for simulation of the polarization observations of debris disks through the Roman Coronagraph.The Stokes parameters for the linear polarization are obtained from the radiative transfer modeling of the debris disk using MCFOST.They are then converted into orthogonal polarization components and propagated through the instrument, incorporating various noise sources and instrumental polarization effects.The observable polarization fraction is estimated from the processed disk images.

GENERATING DISK MODELS USING MCFOST
To illustrate the different modes of operation of Roman Coronagraph, relevant disk observation of two archetypal debris disks systems: the inner ϵ Eridani system will be simulated through the Hybrid-Lyot Coronagraph (HLC) mode, and the HR 4796A system will be simulated through the wide-field Shaped Pupil Coronagraph (SPC) mode.Although the inner disk of ϵ Eridani has never been resolved, we choose to model this system to exemplify the potential of Roman Coronagraph compared to previous coronagraphic instruments.

ϵ Eridani
ϵ Eridani is a star similar to the early Sun, at a distance of 3.2 pc and T*=5100K, M*=0.82M ⊙ and R*= 0.88R ⊙ (Di Folco et al. 2004;Van Leeuwen 2007;Mamajek & Hillenbrand 2008).The outer debris disk around ϵ Eridani has been resolved at infrared and submillimeter wavelengths (Aumann 1985;Greaves et al. 1998;MacGregor et al. 2015;Booth et al. 2017); the inner disk, however, is currently unresolved, and its structural and grain properties have not yet been constrained (Su et al. 2017;Mawet et al. 2018;Wolff et al. 2023).The debris disk has been typically divided into two components: 1) a warm inner disk, which is sometimes hypothesized as an unresolved excess consisting of two narrow belts (1.5-2 AU and 8-20 AU; (Su et al. 2017); 3 AU and 20 AU; (Backman et al. 2009)) and 2) a resolved cold outer disk (55-80 AU (Su et al. 2017) or 90-110 AU (Backman et al. 2009)) imaged with both ALMA (Booth et al. 2017) and other sub-millimeter instruments (Backman et al. 2009;Greaves et al. 2005).We model the inner warm disk with two narrow belts (1.5-2 AU and 8-20 AU) using MCFOST.Dust properties are taken from Su et al. (2017) and shown in Table 1.We use Mie theory (Mie 1908) as the scattering model in MCFOST as it allows for the complete treatment of polarization.As the modeling suggests two separated inner belts, many possible inner structures are possible, which could indicate the presence of companions, and upcoming JWST observations are expected to better constrain the properties of the system; however, the parameters used here provide a physically plausible model to illustrate the sensitivity of the Roman Coronagraph.
The IR excess estimated from the MCFOST modeled spectral energy distribution (SED) is compared with the observed Spitzer-IRS spectrum obtained from Su et al. (2017) and broadband photometry from Backman et al. (2009) as shown in the left panel of Figure 2. The three models shown in Figure 2 use the same parameters from Table 1 except for the grain composition for the inner-most ring, either 100% astrosilicates, 100% olivine, or astrosilicates (50%)+ olivine (50%) to demonstrate that all the three models estimate similar IR excess, SPF.We estimated a higher IR excess than the observed values using 100% amorphous carbon or 100% graphite, while 100% dirty ice shows lower values (not shown here).Among the three dust grain compositions shown in Figure 2, we generated our disk model using the parameters in Table 1.The model IR excesses match reasonably well with the observed IR excess at mid-IR wavelengths.As our disk model contains only the two inner rings, we do not attempt to match the SED beyond 25 microns, which would represent the outer ring.This simplification has negligible impact on the scattered visible light inside 1 ′′ and allows for improved sampling of the circumstellar region of interest.The disk in our simulations is modeled with an inclination i of 34 • and position angle (PA) of 266 • (Booth et al. 2017) for the narrow band filter with a bandpass FWHM of 56.5nm and central wavelength of 575nm.The scattered light and Stokes parameter images are 256 × 256 pixels in size with a pixel scale of 21.84 mas/pixel.Table 1.Parameters used in the MCFOST modelling of ϵ Eridani from Su et al. (2017).As shown in Figure 2, olivines (100%) or astrosilicates(50%)+olivines (50%) can also be used for the inner-most ring.
The model Stokes parameter images (I, Q and U , in units of W/m 2 ) are converted to Jy and further need to be expressed in terms of four orthogonal polarization components for propagation through the instrument.For a perfect instrument (the instrument Mueller matrix is applied later), I 0 , I 90 , I 45 , I 135 can be obtained as, where I corresponds to the total intensity expressed in Jy; these orthogonal polarization components are converted to photons/s using ζ Pup as the reference star with V =2.23 from OS9 ("OS9") sim-ulations, estimating 9.3985×10 8 photons/s at the primary mirror of the telescope.The estimated polarization intensities in photons/s at the primary mirror are shown in the left panel of Figure 3.To measure the impact of instrumental noise (see Section 5 and 6), we will compare the final modeled polarization fraction of the disk with the MCFOST-simulated polarization fraction.To estimate the polarization fraction, one can use p = , but squaring Q and U introduces a systemic bias in low SNR data (Schmid et al. 2006).Therefore, we convert the Stokes Q and U images obtained from MCFOST to Q ϕ and U ϕ , such that the electric field vector direction (in the polarization) is radially oriented with respect to the central star.Following Schmid et al. (2006), this transformation is given by x * , y * corresponds to the pixel location of the central star, and x, y corresponds to all other pixel locations.Q ϕ /I gives the polarization fraction as U ϕ /I becomes negligible; the position angle θ is estimated as 0.5 arctan(U/Q).The Q ϕ , U ϕ , the polarization fraction (p) and position angle (θ) are shown in the right panel of Figure 3.We estimate the maximum polarized intensity of 0.32 mJy/arcsec 2 and the corresponding polarization fraction of 0.37±0.01 in one resolution element (3×3 pixels) in the direction of forward scattering of the disk using the Mie scattering model.In our ϵ Eridani MCFOST modeled disks, we scale the polarized intensity and total intensity to a surface brightness of 0.168 mJy/arcsec 2 per pixel derived from the expected contrast level of 2 ×10 −8 from the non-detection of inner disk of ϵ Eridani in HST observations from Douglas et al. (2024) 3.2.HR 4796A HR 4796A is an A0V star with a bright, inclined debris disk well-studied (distance = 71.9±0.70pc;Van Leeuwen 2007; Prusti et al. 2016) across the optical and NIR wavelengths in polarization using many instruments, Gemini/GPI, VLT/SPHERE, and VLT/NaCo.Hinkley et al. (2009) obtained the first detection of the NIR polarized intensity of the disk at the ansae.The front and back sides of the disk were later resolved in polarized intensity with GPI (Perrin et al. 2015).The improved spatial resolution and smaller IWA identified a brightness asymmetry along the front side of the disk.The data favored an optically thick, geometrically thin model showing a more substantial forward scattering peak at the smallest scattering angles.Using NaCo and SPHERE at the VLT, Milli et al. (2015Milli et al. ( , 2017bMilli et al. ( , 2019) ) detected asymmetry between the northwest and southeast sides of the disk and measured the polarization phase function of the dust in the disk for the first time.These observations were modeled with MCFOST to derive best-fit parameters for grain properties.The observed polarization fraction was found to be 0.4±0.26(40%± 26%) at 90°scattering angle in the optical band (VBB broadband filter), and the observed averaged polarized phase function was compared with both a Henyey-Greenstein phase function model and a Mie theory model with micron-sized dust grains.
Furthermore, multi-wavelength polarization NIR observations presented in Arriaga et al. ( 2020) provided consistent polarization fraction measurements to Hinkley et al. (2009) and Perrin et al. (2015) but showed that single grain composition modeling through Mie (spherical grains) and DHS (distributed hollow spheres) (Jones 1988) theories could not reproduce both the observed polarized and total intensity scattering phase functions (SPF) simultaneously.In a different study, Chen et al. (2020b) used a DHS grain model with a grain composition of silicates (42%), carbon(17%), and metallic iron (37%) to model the VLT/SPHERE H2 total intensity SPF.
At present, the geometrical parameters of the disk derived from all the multi-wavelength observations are largely consistent.In contrast, dust grain models poorly match the observed total intensity phase function and the polarization fraction.This mismatch arises from the shape of the grains and the composition, indicating the importance of using more complicated grain models in addition to improved extraction of the total intensity phase function from the observed data (Tazaki et al. 2019;Arnold et al. 2019).Hence, HR 4796A will be one of the crucial targets to be observed through the Roman Coronagraph in polarization and total intensity, which may better inform the existing dust grain models.
To simulate the polarization observations of HR 4796A through the SPC mode of the Roman Coronagraph, we use one of the best-fit models as an example disk from Milli et al. (2017bMilli et al. ( , 2019) ) shown in Figure 4.The IR excesses estimated for three disk models: a) best-fit SED using Mie theory, b) best-fit SED using DHS theory, and c) best-fit polarization fraction profile (pSPF) are compared with the observed values from Augereau et al. (1999) in the left panel of Figure 4.The polarization fraction obtained from MCFOST is shown in the right panel of Figure 4 and is compared with the VLT/ZIMPOL measurement at 90°scattering angle from (Milli et al. 2019).The best-fit SED models are marginally consistent with the observed polarization fraction.
We use the disk properties from the best-fit SED model given in Table 2 for our simulations.The disk is modeled with an i of 75.8°andPA of 27.7°for the broadband filter with a bandpass FWHM of 96.8 nm and a central wavelength of 825.5 nm.The dimensions and pixel scale of the Stokes images are the same as ϵ Eridani.The Stokes images obtained using MCFOST are converted to orthogonal Best fit SED-DHS (Milli et. al 2017) Best fit SED-Mie (Milli et. al 2017) Best fit pSPF (Milli et. al 2019) Milli et. al (2019) Figure 4. Left: IR excesses for HR 4796A generated with MCFOST for the best-fit models from Milli et al. (2017bMilli et al. ( , 2019)).The observed IR excess is shown from Augereau et al. (1999).Right: The polarization fraction at 825 nm (central wavelength of SPC band) obtained from MCFOST for the three different bestfit models from Milli et al. (2017bMilli et al. ( , 2019)), along with the VLT/ZIMPOL (600-900 nm) measurement of polarization fraction at 90°scattering angle (Milli et al. 2019) polarization components in photons/s using ζ Pup as the PSF reference star (estimating 2.560×10 8 photons/s at the primary mirror of the telescope obtained from "OS11" simulations in the SPC mode of Coronagraph.)The left panel of Figure 5 shows the orthogonal polarization components and the right panel shows Q ϕ , U ϕ , p, and θ estimated Equation 3 and 4. We estimate a peak polarization fraction of 0.92 ± 0.01 and a polarized intensity of 35.25 mJy/arcsec 2 , respectively.The orthogonal polarization components of disk models of ϵ Eridani and HR 4796A shown in the left panel of Figures 3 and 5, have to be convolved with the Point Response Function (PRF) of the Roman Coronagraph instrument which is described in the following section.

GENERATING PRFS FOR THE ROMAN CORONAGRAPH
For each coronagraph mode, a dataset of Point Response Functions (PRFs) is generated using the end-to-end CGI propagation models in roman-phasec-proper (v1.2.5) utilizing PROPER as the back-end propagator.Note that the term Point Response Function is used instead of Point Spread Function as PSF will often imply a linear and shift-invariant instrument response.Due to the influence of apodizers and the focal plane mask, we do not assume a shift-invariant response, so a standard convolution with a PSF cannot be used to generate simulations of disks at the detector.Instead, the dataset of PRFs for a particular mode is interpolated to the array of pixel coordinates of the disk model.This method reduces the image simulation step of a particular disk to a matrixvector multiplication, as explained in Milani & Douglas (2020).Crucially, by utilizing a large set of PRFs that sample the FOV in the radial and angular coordinates, as shown in Figure 6, this method captures the field dependence of the coronagraph PRF within the image simulation.Additional PRFs extending beyond the OWA capture the scattering contributions from sources extending beyond the nominal FOV.
To form the PRF matrix, given a wavelength and source offset, a wavefront is propagated through the optical train to simulate a single monochromatic image.Each HLC PRF is an incoherent sum of seven wavelengths within the band 1 filter centered at 575nm.The SPC PRFs each use five wavelengths within the band 4 filters centered at 825nm.The polarization aberration setting is set to the mean of all polarization states (polaxis=10 within the PROPER models) as we incorporate the polarization effects from the Mueller matrix of the Roman Coronagraph in the final step of simulations.The wavefront of each PRF is normalized to have a total amplitude of 1 at the entrance pupil of the Roman aperture.
The convolved disk images of the orthogonal polarization components I 0 , I 90 , I 45 , and I 135 are shown in Figure 8 for ϵ Eridani (left panel) and HR 4796A (right panel) with the inner working  8, 9.7, 15.2 (140.52, 486.83, 762.87 in mas).The HLC uses a more dense grid of PRFs to reduce numerical artifacts found in the HLC simulations from the interpolation of the PRFs.Right: On-axis PRF corresponding to stellar leakage.The wavefronts of each PRF are normalized at the Roman entrance pupil such that the total sum is 1.
angles (IWAs) and outer working angles (OWAs) of the coronagraphs overlaid in red.The disk and PSF pixel scales are maintained to be consistent during the convolution.

GENERATING RAW EMCCD IMAGES
The Roman Coronagraph will use a back-illuminated electron-multiplying CCD sensor (e2v CCD201-20) consisting of 1024× 1024 pixels of 13 µm in size.It can be operated in low gain (<1000) and high gain (>1000) modes.We use emccd detect (Nemati 2020) to simulate the raw EMCCD images from the convolved disk images.A stack of 50 EMCCD frames is simulated for each orthogonal polarization component with an exposure time of 5s/frame for ϵ Eridani and 1s/frame for HR 4796A, a gain of ≤ 200 incorporating bias of 700e − , dark current of 0.0028 e − /pix/s, and read noise of 100 e − and also incorporate photon noise.Figure 9 shows one of the frames at the EMCCD for all four orthogonal polarization components for ϵ Eridani and HR 4796A, respectively.

INCORPORATING NOISE AND UNCERTAINTY FACTORS FROM OS SIMULATIONS AND DISK PROCESSING
OS simulations are the simulated science images created by the integrated modeling team at NASA-JPL using the most recent version of the observation strategy.They include end-to-end Structural Thermal Optical Performance (STOP) models of the Roman observatory, coronagraph masks, diffraction, wavefront control, detector noise, and jitter.The observing sequence starts with a slew from a servatory settling time before the coronagraphic observations as described in Ygouf et al. (2021).In the most recent OS11 simulations, there are four coronagraphic observation cycles with time between the 2nd and 3rd cycle for the dark hole maintenance.Assuming the dark hole was previously dug and required only minor modifications, each observation cycle begins with observing the reference star (ζ Pup) for 45 minutes, followed by 100 minutes of target star (47 Uma) observations at each of 4 rolls, alternating between -13 °and +13°twice for a total of 400 minutes on target per cycle.The reference star is imaged for 45 minutes again at the end of the cycle.During the four cycles, the reference star is observed six times for a total of 4.5 hours and the target star for 26.67 hours.
The Roman observatory's STOP model is run for a specified timestep to simulate the aberrations and pupil shifts during each observation cycle.Next, the Jitter model produces the RMS jitter over a specified period.Then, the Low order Wavefront Sensing (LOWFS) model is run to generate the Deformable Mirror (DM) correction patterns and is fed to the roman-phasec-proper diffraction model of the observatory.For each timestep, roman-phasec-proper produces complex values of speckle electric fields for four different polarizations (which can be either used individually or two orthogonal polarizations are added for the case of unpolarized source).Finally, these speckle images are propagated through the EMCCD model to incorporate detector noise and uncertainties.Thus, these simulations incorporate all optical aberrations, pointing jitter, DM thermal drifts, polarization aberrations, and EMCCD noise characteristics.These simulations produce data sets with and without noise and also with and without optical model uncertainty factors (MUFs).
In our simulations for the HLC mode, the raw EMCCD images of the disks are rotated to the corresponding roll angles following the steps in the observing sequence, and all the noise components and optical model uncertainty factors (MUFs) from the "OS9" are added.The speckle field images we use in our simulations are [14375,67,67] in dimensions, where each speckle field image is obtained for an exposure time of 5s following the OS time sequence.Similarly, for the SPC mode "OS11" time series speckle field images are added to the raw EMCCD images to incorporate the speckle noise where the speckle field image has dimensions of [1830,181,181] generated for an exposure time of 1s.There are two modes of disk processing to generate the final CCD image as described in the OS simulations: "photon counting" mode with low read noise for gain>1000 and "analog" mode with high read noise for gain<1000.We use the "analog" mode (corresponding to conventional CCD image processing as our disk targets have high SNR) as we use the gain<1000 to simulate all the orthogonal polarization components as shown in Figure 10.In the OS simulations, the target star used is 47 Uma with V = 5.4, and in our simulations, our disk target host stars ϵ Eridani has V = 3.73, and HR 4796A has V = 5.744.We are currently not scaling the speckle fields according to the brightness of our host stars and use the speckle fields generated by the OS simulations as our disks are brighter than the speckle noise.We scaled the I 0 of ϵ Eridani disk a hundred times fainter and processed with the OS simulations to understand the speckle noise level.Figure 13 in Appendix B shows ϵ Eridani processed in the analog mode and the corresponding speckle field noise images.

ESTIMATING POLARIZED INTENSITY AND POLARIZATION FRACTION
The final step in the simulation process is to estimate the Stokes parameters and polarization fraction using the processed disk images of the four orthogonal polarization components.The (Q out , U out ) and total intensity (I out ) are calculated using The instrumental polarization effects, or shifts in Stokes parameters, introduced due to the telescope and instrument optics are represented as a Mueller matrix (Keller 2002).The field-independent Mueller matrix for the Roman Coronagraph is provided by the modeling team5 for wavelengths from 450nm to 950 nm.

ϵ Eridani
We estimated the averaged Mueller matrices for the HLC band and obtained corrected output Stokes parameters, Q cor and U cor as The instrumental polarization, 0.92%, and polarization rotation 0.99×10 −6 obtained from the instrument Mueller matrix are within the expected measurement error and can be easily calibrated.The Lu-Chipman decomposition of the Mueller matrices in Doelman et al. (2023) shows that the negative Q out is obtained due to a mirror and a retarder in the optical path.The Q cor and U cor are converted to Q ϕ and U ϕ using Equations 3 and 4. The "OS9" repository consists of normalized off-axis PRFs for the Roman Coronagraph, which includes losses from the masks but not from reflections, filters, and Quantum Efficiency (QE).We convolved the flux of the reference star ζ Pup with the off-axis PSFs to determine the Zero Point (ZP) magnitude (16.11) to correct for the instrument throughput.
The polarized intensity Q ϕ and total intensity images for ϵ Eridani are shown in the left panel of Figure 11 along with p and θ.We estimate the peak value of Q ϕ as 0.26 mJy/arcsec 2 , I out as 0.78 mJy/arcsec 2 , and a peak polarization fraction of 0.33±0.01 in one resolution element.The peak value of the input polarization fraction shown in the right panel of Figure 3 is 0.37±0.01.Thus, we have successfully recovered the input polarization fraction within the measurement error < 3% after incorporating all the noise sources for a more realistic ϵ Eridani inner disk.

HR 4796A
The averaged Mueller matrices for the SPC band are used to obtain the corrected output Stokes parameters, Q cor , and U cor as, Q cor = −0.005− 0.99Q out + 3.07 × 10 −5 U out (10) Q cor , and U cor are then converted into Q ϕ and U ϕ using Equation 3. The off-axis PRFs from the OS11 repositories are used to estimate the Zero-Point magnitude (18.46), using the reference star ζ Pup for converting ph/pix/sec to their corresponding fluxes.The output polarized intensity and polarization fraction are shown in the right panel of Figure 11.We estimate the peak value of Q ϕ as 32.29 mJy/arcsec 2 , total intensity as 79.03 mJy/arcsec 2 , and a peak polarization fraction of 0.80±0.03 in one resolution element.The retrieved polarization fraction for HR 4796A varies on the order of 0.08-0.10(8-10)% compared to the input polarization fraction.To investigate this discrepancy, We estimated the polarization fraction after each simulation step.The polarization fraction estimated after the convolution showed a difference of 0.08-0.10(8-10)% with the input polarization fraction (shown in Appendix A), which may be a systematic bias causing reduction of the peak brightness of the disk and can be addressed with accurate forward modeling.Thus, it should be an important part of developing the Roman CGI polarization calibration pipeline.However, we have demonstrated that the polarization observations of HR 4796A through Roman-CGI will help accurately measure the polarization fraction and, hence, may help place better constraints on the dust properties.2. We used MCFOST to model two debris disks, ϵ Eridani and HR 4796A, and propagated the orthogonal polarization components through instrument simulation tools.We retrieved the peak polarization intensity and the peak polarization fraction from these simulations.
3. For simulating ϵ Eridani through the HLC mode, using astrosilicates as the dust composition, we recovered the input polarization fraction of 0.33±0.01 at the forward scattering peak after incorporating instrumental polarization and crosstalk.
4. Through the SPC mode, we simulated polarization observations of HR 4796A using the best-fit SED parameters derived from ground-based observations in the optical and NIR.We recovered a peak polarization fraction of 0.80±0.01after incorporating polarization effects from the Roman Coronagraph.
5. We find a difference of ∼0.03-0.10(3-10)% in the output polarization fraction with the input for both of the disks processed using HLC and SPC mode after performing the convolution.This indicates a systematic reduction of the peak brightness of the disk, which must be addressed with accurate forward modeling.The difference between the input and the output polarization fraction is due to the strong PSF smearing effect, which is higher for HR4796A as it's a narrower and sharper ring compared to ϵ Eridani considered here.
6.For the two disks used in our simulations, we obtained sufficiently high SNR with an exposure time of ∼ 250s (5s×50 frames) and hence may not require the target acquisition time of ∼ 26 hours used in the Observing Scenario simulations.Future modeling and simulation efforts are required to derive the optimal exposure times for Roman disk targets.
As a technology demonstration (Kasdin et al. 2020), the coronagraph is no longer bound by scientific requirements.This work, however, validates that the Roman Coronagraph design meets the science requirement developed early in the design process to: "map the linear polarization of a circumstellar debris disk that has a polarization fraction greater or equal to 0.3 with an uncertainty of less than 0.03" (Douglas et al. 2018).This study focused on developing and validating a simulation pipeline for Roman Coronagraph polarimetric observations and demonstration of recovering the polarization fraction from processed disks without considering polarization aberrations (Millar-Blanchaer et al. 2022) and performing "photon counting".The pipeline for the simulated polarization observations of ϵ Eridani and HR 4796A is publicly available (Anche 2023).Future work will incorporate the effects of polarization aberrations from the telescope and the coronagraph, perform photon counting, and compare different post-processing methods (e.g., Karhunen-Lo'eve Image Processing (KLIP) (Soummer et al. 2012) and non-negative matrix factorization (NMF) (Ren et al. 2018)) for disk extraction, and ultimately assess the retrieval of disk geometric and grain properties.The polarization fractions estimated for ϵ Eridani and HR 4796A before and after convolving with the Roman PRFs are shown in Figure 12.The difference of ∼ 3% and ∼ 10% is observed in the case of ϵ Eridani and HR 4796A, respectively.

B. FAINTER DISK INJECTION
In our simulations, we are currently not scaling the speckle fields according to the brightness of the host stars ϵ Eridani and HR 4796A and use the speckle fields generated by the OS simulations as our disks are brighter than the speckle noise.We scaled the orthogonal polarization component I 0 by 100 times fainter and processed with the OS9 simulations to understand the level of speckle-noise fields as shown in Figure 13.The left panel shows the speckle field noise from the OS9 simulations added to the I 0 component.The mid-panel shows the I 0 after the PSF subtraction and the last panel shows the PSF subtracted speckle field image without the I 0 component injected.

Figure 2 .
Figure2.Left: IR excess estimated from MCFOST for ϵ Eridani using astrosilicates (100%), olivines (100%), and astrosilicates (50%)+olivines (50%) for the inner-most ring and H20 dominated dirty ice (100%) for the central ring, compared with observations from(Su et al. 2017) and broadband photometry from(Backman et al. 2009).Right: Corresponding scattering phase function (SPF) at 575 nm (central wavelength of HLC band) obtained from MCFOST for the dust grain compositions.The three dust grain compositions show a reasonable match with the observed IR excess and show similar SPF.

Figure 3 .
Figure 3. Left: The four orthogonal polarization components (I 0 , I 90 , I 45 , I 135 ) in photons/s estimated using the Stokes parameters (Q and U ) and total intensity obtained from MCFOST for ϵ Eridani.Right: Input Stokes Azimuthal components Q ϕ , U ϕ and corresponding polarization fraction Q ϕ /I and θ obtained from MCFOST for ϵ Eridani.A peak polarization fraction of 0.37±0.01 in one resolution element of 3×3 pixels is estimated in the direction of forward scattering of the disk using the Mie scattering model with 100% astrosilicates as the dust grain composition.

Figure 5 .
Figure 5. Left: The four orthogonal polarization components (I 0 , I 90 , I 45 , I 135 ) in photons/s estimated using the Stokes images (Q and U ) and total intensity image obtained from MCFOST for HR 4796A.Right: Input Stokes radial components Q ϕ , U ϕ , corresponding polarization fraction Q ϕ /I, and θ obtained from MCFOST.

Figure 6 .
Figure 6.Left: Distribution of PRFs used for the HLC disk simulations along with the on-axis PSF for the HLC mode.The three concentric rings indicate the IWA, OWA, and maximum radial PRF used for each mode.Respectively, these values, in units of λ/D, are 2.8, 9.7, 15.2 (140.52,486.83, 762.87 in mas).The HLC uses a more dense grid of PRFs to reduce numerical artifacts found in the HLC simulations from the interpolation of the PRFs.Right: On-axis PRF corresponding to stellar leakage.The wavefronts of each PRF are normalized at the Roman entrance pupil such that the total sum is 1.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Left: Distribution of PRFs used for the SPC-WFOV disk simulations along with the on-axis PRF.The three concentric rings indicate the IWA, OWA, and maximum radial PRF.These values, in units of λ/D, are 6, 20, 25.2 (432.06,1440.21,1814.66 in mas), respectively.Right: On-axis PRF corresponding to stellar leakage.Once again, the wavefronts of each PRF are normalized at the Roman entrance pupil such that the total sum is 1.
this work were supported by the WFIRST Science Investigation team prime award #NNG16PJ24 and the Arizona Board of Regents Technology Research Initiative Fund (TRIF).JA is supported by a NASA Space Technology Graduate Research Opportunity.The authors would like to thank Dr. Bruce Macintosh, Dr. John Krist, and Dr Kate L Su for their support and useful discussions.Software: astropy, EMCCD detect, MCFOST, nmf imaging , pymcfost, pysynphot, scipy, pandas PROPER APPENDIX A. DISCREPANCY IN THE POLARIZATION FRACTION AFTER CONVOLUTION

Figure 12 .
Figure 12.The difference in the polarization fraction seen before and after the convolution is shown for ϵ Eridani (top panel) and HR 4796A (bottom panel).

Figure 13 .
Figure 13.Left: The I 0 component of ϵ Eridani is scaled a hundred times fainter and added to the speckle field noise.Center: The I 0 component is processed with the speckle field images.Right: Speckle field image without the I 0 injected.

Table 2 .
Disk properties used in the MCFOST modelling of HR 4796A fromMilli et al. (2017b)for best-fit SED using Mie theory.
. The models for best-fit SED and best-fit polarized intensity scattering phase function (pSPF) predict different IR excess values and polarization, a tension that Roman polarimetry has the potential to resolve.