Wayne - A Simulator for HST WFC3 IR Grism Spectroscopy

Wayne is an algorithm that simulates Hubble Space Telescope (HST) Wide Field Camera 3 (WFC3) grism spectroscopic frames including sources of noise and systematics. It can simulate both staring and spatial scan modes, and observations such as the transit and the eclipse of an exoplanet. Unlike many other instrument simulators, the focus of Wayne is on creating frames with realistic systematics in order to test the effectiveness of different data analysis methods in a variety of different scenarios. This approach is critical for method validation and optimising observing strategies. In this paper we describe the implementation of Wayne for WFC3 in the near-infrared channel with the G102 and G141 grisms. We compare the simulations to real data, obtained for the exoplanet HD 209458 b to verify the accuracy of the simulation. The software is now available as open source at https://github.com/ucl-exoplanets/wayne.


INTRODUCTION
The Hubble Space Telescope (HST) has been invaluable for exoplanetary science. The STIS spectrograph was first used for time-series photometry of HD 209458 b (Brown et al. 2001) and then for the first measurement of an exoplanet atmosphere with transmission spectroscopy (Charbonneau et al. 2002). It has enabled several atomic and ionic detections (e.g. Vidal-Madjar et al. 2003, 2004 as well as scattering properties (e.g. Sing et al. 2011Sing et al. , 2015. The ACS and COS instruments have both made measurements of exoplanet atmospheres (e.g Pont et al. 2008;Linsky et al. 2010). NICMOS was used to search for molecules in an exoplanet atmosphere (e.g. Swain et al. 2008;Tinetti et al. 2010;Crouzet et al. 2012) and has directly imaged exoplanets (e.g. Song et al. 2006). In 2009 the Wide Field Camera 3 (WFC3) instrument was installed and has provided many significant measurements in the near-infrared (IR) channel (e.g. Berta et al. 2012;Swain et al. 2013;Kreidberg et al. 2014b). Whilst WFC3 has been successfully used, as with any instrument, it presents some unique challenges for data analysis. This is particularly true when used in spatial scan mode for exoplanet transit and eclipse spectroscopy, which is the focus of our paper.

Staring Mode
Exoplanet transit spectra were initially observed with HST/WFC3 using the traditional "staring" mode, where the telescope pointing is fixed on the target star (e.g Berta et al. 2012;Swain et al. 2013;Wilkins et al. 2014). The main systematic found in these data is the "ramp" or "hook" effect which causes a sharp rise in flux before quickly plateauing after each HST buffer dump. The other important effect is a gradual reduction in the level of flux throughout the observation. The analysis by Berta et al. (2012) found the amplitude of the "hook" r.varley@ucl.ac.uk to be around 0.4% and the visit long slope to be around 0.05%. These are both much larger than the 10 −4 variations expected in the planetary spectrum. In addition, the initial stage of the analysis includes reductions such as flat-fielding and a non-linearity correction which, despite not being expected to cause issues in the analysis, their end-to-end effect and coupling with other noise sources and systematics has not been fully investigated yet.
Data taken by HST is also non continuous as HST is in a low earth orbit. A typical transit event observed by HST lasts around one to four hours, meaning an observation must span multiple orbits, as a single target is only continuously visible for ∼ 45 minutes. This time includes overheads such as acquiring and re-acquiring the guide star (6 and 5 minutes respectively, Dressel 2016, pg. 209) and buffer dumps. The buffer is a temporary storage space for images on HST which must be transferred (dumped) to the solid state recorder when full. A buffer dump is required after the instrument has taken a certain number of exposures and takes ∼ 6 minutes (Dressel 2016, pg. 212) during which no new exposures are taken. Buffer dumps are generally required when taking many short exposures (usually one or two dumps per orbit for staring data), and lower the overall efficiency of an observation.

Spatial Scan Mode
In 2011 (HST Cycle 19) the spatial scan mode was introduced. This mode slews the telescope during an exposure under the control of either the fine-guidance sensors (FGS, for scans speeds below 5 s −1 ) or gyroscope only control (for scan speeds up to 7 .84 s −1 ). The advantages of this mode over staring include the ability to observe targets that would normally saturate the detector and to achieve higher signal-to-noise ratios (S/N), due to reduced overheads (as each individual exposure can last much longer before saturation). These overheads include the reduced number of detector reads and the reduction, or elimination, of mid-orbit buffer dumps. Many observations have been performed in spatial scan mode (e.g Deming et al. 2013;Knutson et al. 2014a,b;Kreidberg et al. 2014bKreidberg et al. ,a, 2015McCullough et al. 2014;Crouzet et al. 2014;Fraine et al. 2014;Stevenson et al. 2014) however there are some caveats. In addition to the hook present in staring mode, the use of the FGS to control the scan introduces some jitter, causing a periodic "flicker" in the brightness along the scan direction due to FGS feedback (McCullough & MacKenty 2012). Finally, in spatial scan mode there are dispersion variations along the scanning direction (Tsiaras et al. 2016a), shifts in the x and y position of the spectrum in each exposure (Deming et al. 2013) and changes in the effective exposure time when scanning in different directions (the upstream / downstream effect, see McCullough & MacKenty 2012;Knutson et al. 2014a;Tsiaras et al. 2016b).
Since instrumental systematics are of the order of the planetary signal, and as observations typically consist of a single transit, it is hard to verify beyond any reasonable doubt whether the result obtained is in fact the independent planetary signal, not coupled with systematics or residuals from the analysis. This has led to approaches using machine learning algorithms (e.g. Waldmann et al. 2013;Morello et al. 2014) which attempt to extract the signal 'unsupervised'.
In an ideal world we would calibrate our instruments and pipelines with a host of known stable reference stars. In absence of this, we propose here a simulator capable of replicating the measured data, including systematics. Such a tool will help us explore the effects of different instrument systematics and observing techniques on the final scientific result. Moreover, we will be able to understand the coupling between those effects, a study that cannot be carried out based only on real observations, for which the input signal is not known. We can then test the ability of pipelines and reduction methods to recover the planetary signal in different scenarios, which is an important verification step in data analysis. We can, also, validate the stability and reliability of existing data sets, by taking into account the specific configuration that each one has. Finally, such a simulator can be used to explore different possible configurations for future observations and evaluate their effect on the retrieved planetary spectrum.
To this end, we describe here the implementation of Wayne, a simulator for WFC3 IR spectroscopy in both the G102 and G141 grisms. We implemented many systematics and noise sources present in real data and Wayne is able to simulate both ordinary stellar spectra as well as exoplanet transit spectroscopy, where the planetary and stellar spectra are combined together in a transit event. We demonstrate the accuracy of our simulation by replicating the real observation of HD 209458 b from HST proposal 12181 and process it using the same pipeline described in Tsiaras et al. (2016a). We compare the real observation to our simulation providing additional verification of the analysis for that observation.

Other Simulators
Comparisons are obviously drawn to aXeSim 1 , the STScI grism simulator package for WFC3. aXeSim was designed for general astrophysics whereas Wayne has been conceived for transit and eclipse spectroscopy. The limitations of aXeSim, which are accounted for in Wayne, include spatial scan mode, scan speed variations, dispersion variations, x and y positional shifts, trends such as the "hook" and long term ramp, cosmic rays and the ability to simulate a transit with time. However, Wayne is currently only able to simulate a single spectrum at a time whereas aXeSim can simulate a whole field at once. We also do not attempt to reimplement the functionality of the Hubble Space Telescope Exposure Time Calculator 2 , which is able to approximate the signal to noise of a target. We instead rely on the user to scale the input flux of the star to the correct level using the exposure calculator or otherwise. The software is now available as open source at https://github.com/ucl-exoplanets/wayne.

Wayne
The simulator imitates real data taken by WFC3 IR grisms through four key steps: • Preparing the visit -Calculation of the exposure time of each frame including any gaps during orbits and buffer dumps.
• Preparing the spectrum -Scaling the flux to the instrument, applying the grism sensitivity, any re-binning and, in the case of transit spectroscopy, combining the stellar and planetary signals.
• Converting flux into pixels -Calculating where on the detector a specific wavelength falls, simulating the point-spread function and integrating it into pixels.
• Adding any noise and implementing reductions -Applying the flat field, scan speed variations, detector gain, read noise etc.
These steps are often overlapping as reductions and systematics need to be applied at different stages of the simulation. A simplified overview of the process used to create a simulation is shown in Figure 1.
There are many sources of noise and trends caused by the telescope and its instruments. Some of these, such as the dark current, flat-field and quantum efficiency have been constrained very well through calibration programs and are used in reduction pipelines like aXe (Kümmel et al. 2009). Others, like the variation in scan speed and the hook trend are not yet well understood or calibrated. For the first set we used the state of the art from the literature. For the latter, we conducted our own investigations on real data sets to determine approximate models for these features, and allow their parameters in the simulations to be varied by the user. We plan to investigate these behaviours in the future. The WFC3 detector is read multiple times per exposure in what is known as a non-destructive read or "sample up the ramp". We refer to these as "reads". In order to generate a spatial scan exposure we must combine many staring mode frames together at intervals in time (and therefore position). We use the term "subsample" to refer to these to avoid confusion. The exposure time of a subsample is set in the configuration file of the simulator and is generally 5-30 ms. The exact, i.e computationally optimised, value of this parameter depends on the scan speed.
In this section we mention certain files such as the 'flatfield cube' which contains the corrections we are using. These are mostly obtained from STScI and we provide the location for each file in the appendix B.
In this paper we use the transit spectroscopy mode while spatial scanning as an example of Wayne. The G102 and G141 grism implementations are very similar, differing only in the calibration files used. As such, we focus on the G141 grism as it is more commonly used in exoplanet transit spectroscopy.
2.1. Planning a Visit HST is located in a low earth orbit, completing an orbit of Earth in around 95 minutes. Single targets are typically visible for 50 minutes of each orbit but can be longer depending on their location with respect to HST orbit. Visits comprise of multiple orbits, typically limited to 5 or 6 due to the South Atlantic Anomaly (Gonzaga 2015, pg. 17). The first orbit in a visit exhibits the most severe trends due to spacecraft settling and is therefore normally left out of the analysis. For this reason, we do not attempt to recreate the first orbit effects in our simulations.
WFC3 IR detector has several different observing modes which set the exposure time and the size of the detector array that is read. These parameters are defined by three variables, the subarray mode (SUB-ARRAY =1024, 512, 256, 128, 64), the sample sequence (SAMPSEQ=RAPID, SPARS5/10/25/50/100/200, STEP25/50/100/200/400), which describes how the non-destructive samples up the ramp are spaced and the number of samples (NSAMP =1 to 15). The exposure time can be obtained for a particular combination from Rose (2015, section 13.3.6). Wayne is able to simulate all permitted modes, we are only limited by the available configuration files and exposure time information given.
Wayne's visit planner is a function that takes all the information about the visit and calculates the exposure start times, orbit ends and buffer dumps. Specifically it takes input of NSAMP, SUBARRAY, SAMPSEQ, the number of orbits and whether the telescope is scanning. It then outputs the start time of each exposure in each orbit. This involves taking into account target visibility, guide star acquisition or re-acquisition, exposure time, read time, spacecraft manoeuvres and buffer dumps from information given in Dressel (2016, chapter 10). The produced time sequences do not match perfectly real data or the observing schedules from the Hubble Space Telescope Astronomer's Proposal Tool 3 , due to the approximate values given in the handbook. When recreating a real observation, the exposure start times from the real data can be given allowing for a more precise simulation.

Non-dispersed Calibration Image
For wavelength calibration in the analysis of WFC3 spectroscopic data a non-dispersed (direct) image of the star is taken at the beginning of an exposure. This is used to calculate the x and y reference positions. We output a simple observation for this file consisting of a 2D Gaussian centred on the reference position to simulate the star. This allows the simulation to be analysed by existing pipelines that require this frame.

Transit Spectroscopy: Combining Signals
When used in transit spectroscopy mode, the simulator must combine the planetary and stellar signals together. To do this, we need to simulate a transit event. When a planet transits a star it blocks a fraction of the stellar light we receive. As the planet moves across the stellar disc, we observe a light-curve. In order to simulate a transit we need a model of the stellar spectrum (the flux as a function of wavelengths), the planet spectrum (the transit depth as a function of wavelength) and then combine them at time t and wavelength λ. A transit light curve is required for each spectral bin. They are calculated based on the above input signals and the orbital parameters of the planet, using our PyLightcurve 4 package (in prep.) which implements a numerical light curve model. The advantage of this implementation over other models (Mandel & Agol 2002;Giménez 2006;Pál 2008, e.g.) is flexibility in the use of the non-linear limb darkening law (Claret 2000), support for eccentric orbits, and it is written in pure Python (ensuring compatibility with the simulator). An example of the light-curve model for HD 209458 b (described later) is shown in Figure 2.

Flux to Pixels
The wavelength-dependent flux from the star passes through the grism and is distributed to the detector pixels. This involves scaling the flux by the grism throughput and the quantum efficiency of the detector. Then, by using the wavelength calibration coefficients of the grism we find the centre of the point-spread function (PSF) on the detector and distribute the flux to pixels based on the PSF morphology at that wavelength. We describe each of these steps in more detail below.

Scaling Flux
For the simulation we assume the input stellar flux is pre-scaled to the correct value for the WFC3 IR channel. We provide an option to add a scaling factor which can be used to manually adjust the flux to this value. We do not, at present, attempt to implement the functionality of the Hubble Space Telescope Exposure Time Calculator.

Grism sensitivity
The sensitivity of the grism was calibrated by Kuntschner et al. (2011a) for each order of the spectrum. We only consider the first order spectrum in this version of the simulator, as it is the main scientific spectrum.   The sensitivity is given in units of cm 2 e − erg and quantifies the conversion of flux to electrons, accounting for both the throughput of the grism and the quantum efficiency of the detector. As shown in Figure 3, the sensitivity of the grism doubles between the G141 wavelength limits 1.1 to 1.7 µm. We use a linear interpolation to scale the sensitivity to the spectral bins.
The quantum efficiency (QE) of the WFC3 IR detector is less significant than the grism throughput but is still an important effect. The QE was calibrated on-orbit (Baggett et al. 2010) and is shown in Figure 4. The QE is largely flat for the G141 grism but has a significant uptrend for the G102 grism. We do not apply a separate QE step in our simulation as it is contained in the grism sensitivity.

The field-dependent structure of the spectrum
When a grism is used, the source is spread out into a spectrum. The line through this spectrum is known as the trace line, while the equation giving the wavelength as a function of distance on the trace line is know as the wavelength solution. According to Kuntschner et al. (2009a,b) both the trace of the spectrum and the wavelength solution depend on the reference or source position on the detector (field-dependent). This is the central location of the observed star on the detector without the grism and is often obtained from a single non-dispersed (direct) image taken at the beginning of the visit.  Using the field-dependent wavelength solution given by Kuntschner et al. (2009a,b) we can map each pixel to a wavelength. An example trace for a given source position is given in Figure 5. According Where x and y are the detector positions of the trace and the coefficients a n are defined by Kuntschner et al. (2009a,b) for both grisms ( Table 6).
The wavelength solution is given by where λ is the wavelength, d is the distance between the reference and the required point along the spectrum trace and b n are the coefficients defined by Kuntschner et al. (2009a,b) for for both grisms (Table 7).

PSF
The point-spread function (PSF) for the WFC3 IR channel can be modelled as the linear combination of two 2D Normal distributions (N) with different standard deviations (σ): aN 1 (σ 1 ) + (1 − a)N 2 (σ 2 ). The values for the three different parameters are modelled based on the PSF ensquared energy fraction as a function of aperture size and wavelength (Dressel 2016, pg. 144).
The PSF parameters are shown in Figure 6 as 3 rd order polynomial functions of wavelength, a dependency that significantly affects the final spectrum (for the coefficients used see Table 8 (pixels) Figure 6. The PSF parameters with wavelength for the WFC3 IR channel. can be found in Table 1. The difference is due to the diagonal components of the HST PSF, which we do not include. However, the size of the PSF wings is comparable with the real data. The final step is to distribute on the detector the total number of electrons for each wavelength bin, given the centre and the shape of the PSF at this wavelength. We associate to every electron two random numbers drawn from one of the two normal distributions that construct the PSF -a fraction of 1/a of all the electrons belong to N 1 and a fraction of 1/(1-a) belong to N 2 . These random numbers are then added to the centre of the PSF to give the x and y pixel coordinates of all the electrons. To optimize this process, the random number generation is implemented in C, using the Box-Muler method, and is accessed by Python using the Cython compiler (http: //cython.org/). At this point, we have described all the parts needed to form a simple staring mode frame. Using the parameters described in section 3 we have generated a staring mode frame shown in Figure 7.

Noise and data reductions
There are many different reductions needed for processing WFC3 data which are fairly well calibrated from on-orbit calibration studies. Noise sources have been quantified in a similar way. We describe here how we implement each of these steps in the simulation. Note that in our figures we show calibrations for the 256 subarray as this is one of the most commonly used for exoplanet transmission spectroscopy with WFC3.

Photon Noise
Photon noise is added by converting the stellar flux rate for each sub-sample into counts using a Poisson distribution before distributing the flux into pixels.

Flat field
For the flat field we use the flat-field cube described in Kuntschner et al. (2011b) (Figure 8). The flat field for grism spectroscopy is a 3 rd order polynomial with the coefficients given in the flat field cube. The flat is calculated on a per subsample basis by determining the wavelength of every pixel at the subsample reference position.

Sky Background
The sky background is generated by multiplying the master sky image (Kümmel et al. 2011) (Figure 9), by a factor depending on the target and then sampling each pixel from a Poisson distribution to add noise. We acknowledge the recent source-dependent master sky images by Brammer et al. (2015) and plan to implement them in a future version.

Cosmic Rays (CR)
Cosmic rays impact the detector at a rate of 11 ± 1.5 s −1 across the full frame (Dressel 2016, pg 153). We add cosmic rays to each sample by scaling the rate by the exposure time and subbarray size and sample a Poisson distribution with this rate to give the number of impacts. We then randomly choose a pixel for the impact to occur and the number of counts between 10-35 kDN to add to that pixel.

Gain variations
The gain is the e − /DN conversion of the amplifier used to read a pixel. This is composed of a constant factor 2.35 and gain variations per quadrant which are specified in the gain calibration file u4m1335mi pfl.fits. We therefore divide the simulation counts by 2.35 and then multiply by the gain variations.  Figure 10. WFC3 IR superdark u4819493i drk.fits for a 256 subarray at NSAMP 5 in SPARS10 mode covering 29.6 s. The frame is scaled between -10 and 5 electrons per pixel. The mean is -1.65, the median is -3.37, the minimum is -14.3, and the maximum is 14000 electrons per pixel. This includes any bad pixels (cold and hot). lude et al. 2014). Despite this, for completeness we add dark current, using the superdarks (Dulude et al. 2014), created by averaging many dark exposures from cycles 17, 18, 19 and 20. They are specific to the subarray and readout mode used, as the dark current scales nonlinearly with time. An example superdark is shown in Figure 10. Dark current is added to each sample at the end of the simulation on a per read basis.

Non-linearity
The WFC3 IR detector reaches 5% non-linearity at 70,000 e − or ∼ 28,000 DN with a maximum well depth of 40,000 DN (Hilbert 2008).
The form of the non-linearity is described by Hilbert (2008) on a per quadrant basis. We note the more precise non-linearity correction per pixel by Hilbert (2014) but the coefficients to implement this correction were not readily available. The per quadrant version is described by a 3 rd order polynomial with the coefficients c given in the calibration file. The conversion from the non-linear F to the linear case F c is then: For our simulations, we want the reverse, i.e. to convert linear counts to non-linear. To do this, we use the Newton-Raphson numerical method to find the root of the above equation that is closest to the original counts.

Count level clipping
As the non-linearity correction only applies to counts up to the 5% non-linearity limit we clip the maximum values of all pixels to this level (70,000 e − ) and consider them saturated.  Figure 11. The initial bias frame for the 256 subarray mode used for the zero-read. The frame is included mostly for consistency with real data as it is removed in zero-read subtraction -the first stage of analysis. The scale has been clipped between 5,000 and 19,000 DN.

Bias and bias drifts
When the detector is reset there is an initial bias level of around 11,000 DN with small differences per quadrant. This bias level gives the look of a raw WFC3 frame before subtracting the zero-read and is nearly fully corrected by subtracting the zero-read (Deustua 2016, pg. 126). For consistency with real data we include this effect by creating an initial bias file from the zero-read of two frames of real data obtained from the STScI MAST Archive 5 . Specifically we use the upper half of the file icdp02b0q raw.fits and lower half of the icdp02b1q raw.fits file (to obtain areas not contaminated by the spectrum). The resulting file is shown in Figure  11.
The bias level drifts during an exposure. The WFC3 IR detector allows correcting for this effect by having a five pixel outer border around the 1014 × 1014 array that is not sensitive to incoming light. These are known as reference pixels and are a proxy to the drifting bias level. By averaging them and subtracting the result from each read we can correct the bias drift. We do not attempt to implement bias drifts in the current simulation, instead we set all the reference pixels to 0 before applying the read noise.

Read Noise
Read noise for the WFC3 IR channel is Gaussian with a standard deviation between 19.8 to 21.9 e − (Hilbert & McCullough 2009) after subtracting the zero-read. For a single read, we therefore sample a normal distribution for each pixel with standard deviation of 20/ √ 2 = 14.1 e − (6 DN) and apply it non-cumulatively to each sample.

Spatial Scanning
Spatial scanning mode is created by combining multiple staring modes frames together. We generate subsamples at the subsample rate, typically 5-30 ms depending on the scanning speed. The planetary and stellar spectra are combined per subsample so the transit is progressing throughout the scan. Trends such as the scan speed variations are applied on a subsample basis.

Reading the array
The WFC3 IR detector is read separately in quadrants with each quadrant read towards the centre of the frame. This process takes 0.278 s in total for the 256 subarray meaning that the flux content of the different pixels is not recorded at the same time. For all pixel exposure times to be equal, the zero-read must be subtracted.
The read time and direction has an impact when scanning due to the up-stream / down-stream effect identified by McCullough & MacKenty (2012). This causes a change in exposure time depending on whether the scan is with or against the read direction. We do not simulate this effect at present but plan to include it in a future version as it has a particular effect when a visit is composed of scans in both directions (see Knutson et al. 2014a;Tsiaras et al. 2016b)).

Samples up the ramp
During an exposure, the detector is read several times, in what are known as 'samples up-the-ramp' which we refer to as reads. We generate each read separately to recreate this data structure. The zero-read is currently implemented as the initial bias described in Section 2.4.9 with read noise.

Output to FITS format
The simulations output to FITS files in the same format as HST WFC3 data and contain most of the same header information. In addition we include several new header keys giving details about the simulation such as the simulator version, x-shift amount and the planet and star parameters used. A full list is given in table A.

Systematics
Data obtained with WFC3 include several systematics that can be measured per exposure over an entire visit. These include the shift in the starting x and y position in each frame (x-shifts and y-shifts), scan speed variations and the 'hook' or 'ramp' effect. the total shift is commonly only ∼ 0.1 pixels (Tsiaras et al. 2016a). However, there are examples where the shift is higher, i.e. a total of ∼ 1.3 pixels in Tsiaras et al. (2016b) where y-shifts can not be ignored. In addition to the position shifts from one frame to an other jitter noise is also present in HST observations. Jitter occurs due to the imperfect pointing of the telescope during an observation, and results in an effective increase in the size of the PSF. Jitter noise, across both axis, is included as uncertainty in the reference position (x ref , y ref ) of the spectrum. To match real data sets, the input value of the jitter noise can be derived from the available jit files. For HST, the jitter is typically between 0.003 and 0.005 (0.022 to 0.037 pixels). For the HD 209458 b data set used in section 3 the deviation from a straight line scanning trajectory, which is assumed in Wayne, is ±0.0013 (0.01 pixels). However, these measurements represent time intervals of 3 seconds, corresponding to 0.0013 √ 3 = 0.0023 (0.017 pixels) for intervals of 1 second. In section 3 we used a jitter of 0.02 pixels.

Flux distribution / positioning
Detector effects included flux-to-electron conversion zero-read random sampling of the PSF read noise wavelength-dependent PSF dark current field-dependent spectrum trace gain conversion field-dependent wavelength solution non-linearity spatial scanning wavelength-dependent flat-field scan speed variations field-dependent sky background position shifts cosmic rays jitter noise not included "blobs" image persistance satellite trails cross talk inter-pixel capacitance charge diffusion non-ideal A/D converters (Section 2.4.9) "snowballs" up-stream/down-stream effect (Section 2.5.2)

Scan speed variations (SSVs)
If the FGS is used for controlling the scan there is a flicker which appears as a wave-like variation in flux in the direction of the scan, due to FGS feedback (McCullough & MacKenty 2012).
We looked at SSVs in the HD 209458 b data set used in section 3. After applying all reductions to the frames we choose a single orbit (orbit 3) and summed the rows of the spectrum to give the total flux per row in the scan direction for each of the 20 frames. Visually these variations appear to be modulated sine waves (see Figure 12 top). We confirmed this behaviour by performing Principle Component Analysis on these 20 summed rows and recovered the first 2 principal components which appear to be 1) the general increase in flux across the rows 2) a sine wave with varying amplitude (see Figure 12 middle). We fit the second principle component to give an approximate sine model of period=0.7 seconds, standard deviation=1.5 and phase=0 (see Figure 12 bottom). These values are configurable in the simulation to adapt to other data sets.
The sinusoidal model is applied to the simulation by changing the exposure time of each subsample to that defined by the sine model with mean equal to the subsample time. We note that the real scan speed variations are reasonably complex, appear to modulate and have occasional large 'blips' in the data consisting of a large peak and trough in the flux level. We only give an initial approximation of the effect for this data set here as an in depth analysis is beyond the scope of this paper.

Hook and long term ramp
There are two more trends often seen in real data which are often combined into a single reduction step. This is the hook effect which appears to reset after each buffer dump and a long term ramp which is a gradual decrease in the level of flux with time across a visit.
The 'hook' was first described by Berta et al. (2012) and implemented as a exponential fit by Kreidberg et al. (2014b). We adopt the form of Tsiaras et al. (2016a) (Equation 8) and apply it as a scaling factor to the flux for each subsample.
Where t is time, t v the time when the visit starts, t o the time when the orbit in which the frame belongs starts, r a the slope of the linear long-term 'ramp' and (r b1 , r b2 ) the coefficients of the exponential short-term 'ramp'. The values of these coefficients can be set in the configuration to match a similar observation. Table 2 summarizes the effects related to the spatial scanning technique and the WFC3 IR detector, and here we discuss the effects that are not currently implemented in Wayne.

Summary of simulated effects
"Bolobs" -The "blobs" are areas of 10 to 15 pixels that absorb up to 15% of the incoming light at their centers. The behaviour of the "blobs" appears to be datedependent, making their implementation in Wayne more difficult than the offer effects (Dressel 2016, pg. 160).
Satellite trails -Such trails can appear when satellites cross the observed field of view. Due to their brightness, satellites saturate a number of pixels (approximately 10 to 15), depending on the length of the trail. Saturation can also affect consecutive observations because of image persistence. In the case of exoplanetary observations, the exposure times are short and the satellite trails are more rare.
Image persistance -Due to image persistence, an exposure of the detector is causing ghost images or afterglows in the consecutive ones (Dressel 2016, pg. 154). This effect is, possibly, the cause of the "hook" seen in each HST orbit. Implementing this effect requires the coupling between different exposures, which is, currently, not available in Wayne. Crosstalk -Is the effect when a bright source in one quadrant is causing electronic ghosting in another quadrant, coupled to the first one. In the WFC3/IR channel, the two quadrants on the left side of the detector are coupled, and the the two quadrants on the right side are also coupled. However, this effect is below the background noise for unsaturated sources (Dressel 2016, pg. 66).

Inter-pixel capacitance (IPC) and charge diffusion (CD) -
Both effects cause a fraction of the charge collected by an individual pixel to "leak" into its neighboring pixels and therefore the recorded image to appear smoother. We plan to include both effects in the future, by using appropriate smoothing kernels. IPC is a form of crosstalk and CD is caused during charge transfer through the detector, hence, only the later is wavelength-dependent (Hilbert & McCullough 2011).
"Snowballs" -"Snowballs" have the same behavior as the cosmic rays, but they affect about 10 pixels and they contain between 200,000 and 500,000 e − . Because they appear only about once per hour over the entire WFC3 detector, this effect is not common in spatially scanned data sets (Durbin et al. 2015).

HD 209458 b: A CASE STUDY
As an application, we simulate the scanningmode spectroscopic observations of the hot Jupiter HD 209458 b (ID: 12181, PI: Drake Deming). The simulated data set consists of 83 images and each image of 5 up-the-ramp samples with a size of 266 × 266 pixels in the SPARS10 mode. As a result, the total exposure time, maximum signal level and total scanning length are very similar to the real frames (22.32 seconds, 4.8 × 10 4 electrons per pixel and 170 pixels, respectively). Our simulation also includes one non-dispersed (direct) image of the target, at the beginning of the observation, for calibration reasons, similar to the real data set. A raw frame from the simulator and real data is shown in Figure 3.
For this simulation we use a PHOENIX model (Allard et al. 2003;Baraffe et al. 2015) for the stellar spectrum for a star as similar as   Figure 13. Comparison of a real raw frame (ibh726meq raw.fits) of HD 209458 from proposal 12181 (upper) and the Wayne simulation (lower). Both represent the flux difference between the last and the first reads and are cropped to the same limits to highlight the spectrum. The main visual difference is the non-uniform scan speed variations in the real frame. Simulated cosmic rays are random and will differ from the real frame.
The total execution time for a single spatially scanned image is of the order of 37 seconds, but it is highly depended on the hardware. Among the different stages of the simulation, 16% of the time is spent on calculating the flux as a function of wavelength, 68% on randomly distributing the flux into the pixels, 9% on applying the wavelength-dependent flat-field and 7% on all the other processes.

Reduction-Calibration-Extraction
We analyse the simulated data set using our specialised pipeline for reduction, calibration and extraction of scanning-mode spectroscopic data (Tsiaras et al. 2016a) for a one-to-one comparison with the real data set.
At first, we apply all the reduction steps in the same way as implemented for a real data set, including: bias drifts correction, zero-read subtraction, non linearity correction, dark current subtraction, gain variations calibration, sky background subtraction and bad pixels/cosmic rays correction.
Then, the frames are calibrated for position shifts. For the (horizontal) x-shifts, we are making use of the normalised sum across the columns of the frames, while for (vertical) y-shifts we determine the position of the spectrum on the first non-destructive read. The input values   Figure 16. 1D spectra extracted from a real and a simulated spatially scanned frame, compared with the sensitivity curve of the G141 grism. Although the stellar spectrum used is the closest to HD 209458, we can see that the real spectrum is smoother, probably due to a lower number of absorption lines.
that we use match the shifts presented in the real data set and, as it can be seen in Figure 14, we are able to retrieve them very accurately. This result ensures that all the effects caused by the position-dependent systematics are properly simulated by Wayne.
To extract the 1D spectra from the spatially scanned spectra we first calculate the physical position of the star on the full-array detector. We then shift the position of the photons during the scan, and extract the flux for each wavelength bin from an aperture of quadrangular shape (Tsiaras et al. 2016a). In Figure 15 we can see the wavelength-dependent photon positions as calculated for a real and a simulated frame. From the similar structure we can conclude that the field dependent spectrum trace and wavelength solution are correctly implemented during the scanning process.
We, can therefore compare the 1D spectra, as extracted from a real and a simulated frame ( Figure 16). The two spectra have the same shape, consistent with each other and the sensitivity curve of the G141 grism. However, the real spectrum is smoother, probably due to a lower number of absorption lines.

Fitting the white and spectral light-curves
The last step is to fit the white-light curve for the midtransit point and the 'ramp' coefficients and then use these values to fit the R p /R * for the different wavelength channels. We follow the same two approaches as we did for the real data set, i.e fitting the instrumental systematics function, Equation 8, and the transit model both separately and simultaneously, and keeping inclination, i, and a/R * ratio fixed. The results we obtain (Table  3) are consistent with each other and in good agreement with the inputs. Also the magnitude of the errors are as expected for the signal-to-noise ratio of this particular target observed by WFC3.
We then fit the spectral light-curves leaving as free parameters only the normalisation factor and the R p /R * ratio. As explained in Tsiaras et al. (2016a), the spectral light-curves include an additional kind of systemat- ics which originates from the low resolution of the spectrum (under-sampling) coupled with the horizontal shifts (Deming et al. 2013). Since the simulated data set includes these shifts we have to take these systematics into account. The final results can be seen in Figure 18. The spectrum extracted from the simulated data has a mean uncertainty of 38 ppm, slightly better than the real one (40 ppm). The rms of the residuals between the input and the output spectra is also 38 ppm and it is consistent throughout 10 different simulations. Based on these results we can conclude that simulations created by Wayne can reproduce existing data sets in a realistic way.

CONCLUSION
We implemented the Wayne algorithm as a Python package capable of simulating Hubble Space Telescope (HST) Wide Field Camera 3 (WFC3) grism spectroscopy, including sources of noise and systematics. From all our diagnostics we conclude that given the configuration of a real data set, Wayne generates a simulated data set with the same behaviour and noise level as the real case. Therefore, it is a powerful tool with which we can validate the stability of existing data sets with certain configuration but also predict the behaviour of future observations in order to decide the best observing strategy for a particular target.  Wayne can then be used to test how different analysis methods work in different scenarios in order to validate data analysis pipelines. We demonstrate in this initial investigation that the method of Tsiaras et al. (2016a) can accurately recover the independent planet signal for that data set.