Dynamical Surface Imaging of lambda Andromedae

We present temperature maps of RS CVn star lambda Andromedae, reconstructed from interferometric data acquired in 2010 and 2011 by the MIRC instrument at the Center for High Angular Resolution Astronomy Array. To constrain the stellar parameters required for this imaging task, we first modeled the star using our GPU-accelerated code SIMTOI. The stellar surface was then imaged using our open source interferometric imaging code ROTIR, in the process further refining the estimation of stellar parameters. We report that the measured angular diameter is 2.742 +/- 0.010 mas with a limb-darkening coefficient of 0.231 +/- 0.024. While our images are consistent with those of prior works, we provide updated physical parameters for lambda Andromedae (R_star = 7.78 +/- 0.05 R_odot, M_star = 1.24 +/- 0.72 M_odot, log L/L_odot = 1.46 +/- 0.04).


INTRODUCTION
Observing stellar surfaces provides insight to the physics within stellar interiors. We know that stars ranging from pre-main sequence to giants exhibit magnetic spot activity on their surfaces (Strassmeier 2009). Since the advent of space missions, such as CoRoT (Baglin et al. 2006a,b) and Kepler Koch et al. 2010), many more stars have been observed to exhibit magnetic activity (Frasca et al. 2011;Fröhlich et al. 2012;Roettenbacher et al. 2013Roettenbacher et al. , 2016aNielsen et al. 2019;Santos et al. 2019). These stellar features constitute major sources of uncertainty trying to calculate accurate stellar physical parameters (e.g., T eff and R ; Somers & Pinsonneault 2015). Starspots have other astrophysical significance tying them to accurately determining exoplanetary parameters. Any uncertainties found in the host star's physical parameters are amplified to any of their planetary parameters, as deriving exoplanetary parameters are dependent on the parent star.
RS Canum Venaticorum (RS CVn) variables are known to present large magnetic starspots (Hall 1976;Kővári et al. 2015;Roettenbacher et al. 2016bRoettenbacher et al. , 2017. These variables are often found in a binary system and the pair often consists of an evolved giant primary with the secondary being a smaller main-sequence companion. Magnetic spots in these systems are often easier to observe because of their relative size to the star, thus making RS CVn variables ideal targets. There are three main techniques routinely employed to image these systems: light-curve inversion, Doppler imaging, and interferometric imaging. Photometric monitoring of these systems provides straightforward evidence for stellar spots, as shown in many systems observed by the Kepler spacecraft (e.g., Frasca et al. 2011;Fröhlich et al. 2012;Roettenbacher et al. 2013Roettenbacher et al. , 2016a. The inverse problem of imaging the stellar surface from photometry is called light-curve inversion (Wild 1989;Roettenbacher et al. 2013). A main drawback of broadband light-curve inversion is that photometry only provides relative information about the latitude of starspots (Harmon & Crews 2000) and relies on a prior knowledge of the stellar limb-darkening. Light-curve inversion from multi-band photometry alleviates the latitude ambiguities, hence resulting in more accurate solutions (Harmon & Crews 2000).
Doppler imaging (Goncharskii et al. 1977; Rice et al. 1981) is the class of inverse methods for imaging stellar surfaces from spectroscopic data. This technique uses perturbations of absorption features on a star to better estimate the spot's latitude and longitude. However, there are still uncertainties in determining spot location for stars near edge-on rotation. High-resolution spectra are needed in Doppler imaging to distinguish the features due to the starspots in the absorption lines and to be able to accurately detect their locations. High rotational velocities rotationally broaden absorption lines and are required to ensure that the spectroscopic impact of a spot moving across the surface is shorter than the spot's evolution timescale. Piskunov & Wehlau (1990) determined lower bounds enabling Doppler imaging to be from 6 km/s to 15 km/s, which corresponds to spectrograph resolving powers of at least 20,000 to 50,000.
Contrary to Doppler imaging or light-curve inversion, interferometry provides unambiguous evidence that a spot is being shown without any assumptions on latitude. Interferometric modeling allows the determination of angular parameters, such as the inclination or position angle of a spotted star. However, interferometric observations can only be managed on a limited number of targets (i.e., relatively bright targets) when compared to photometric and spectroscopic targets, and furthermore only targets of sufficient angular size can be resolved from Earth. It was only in 2007 that interferometric synthesis imaging became possible (Monnier et al. 2007) thanks to longer baselines and the combination of light from four (and now up to six) different telescopes.
λ Andromedae (HD 222107; hereafter λ And) is a bright G8III-IV RS CVn variable (V = 3.82, H = 1.40) with spots, and is included in the third edition of the Catalog of Chromospherically Active Binary Stars (Eker et al. 2008). It is a single-lined spectroscopic binary system with a rotation period of 54.07 days for the primary  and has a companion in asynchronous rotation. Walker (1944) found a nearly circular orbit for the system with an eccentricity of e = 0.084 ± 0.014 and an orbital period of 20.5212 ± 0.0003 days. The most recent estimate of the effective temperature and mass for the primary star of λ And is 4800±100K and 1.3 +1.0 −0.6 M (Drake et al. 2011) with its companion most likely being a low mass main sequence star or a massive brown dwarf based on the mass ratio calculation of q = 0.12 +0.07 −0.04 (Donati et al. 1995). Parks et al. (2021) was the first to do 2D snapshot interferometric imaging of λ And using data obtained with CHARA. Their study estimated the angular diameter for the primary of λ And to be 2.759 ± 0.050 mas, corresponding to a physical radius of 7.831 +0.067 −0.065 R given the Hipparcos distance of 26.41 ± 0.15 pc (van Leeuwen 2007).
In this paper, we describe the process we followed to obtain a temperature map of the surface of λ And. In section 2, we present data acquisition and reduction. In section 3, we describe how we used the interferometric modeling code SIMTOI to obtain initial guesses of stellar parameters. We introduce the ROTIR imaging code in section 4, then its application to the imaging of λ And in section 5. The imaging results are compared with previous works in section 6. We go on to discuss prospects beyond solid body rotation in section 7 and the search for the companion of λ And in section 8. Finally, we discuss our conclusions and future work in section 9.

OBSERVATIONS
We reuse the 2010 and 2011 data from Parks et al. (2021), shown in Table 1 and calibrators in Table 2 used for each respective year, for our analysis. These data were obtained using the CHARA Array (ten Brummelaar et al. 2005) with the Michigan Infra-Red Combiner (MIRC; Monnier et al. 2004) in H-band with the median wavelength of 1.65 µm. The observations were done in prism mode (R = 50) which contain eight spectral channels. The data taken in 2010 were taken with a combination of four out of six telescopes which provide six visibilities, three independent bispectrum amplitudes (triple amplitudes), and three independent bispectrum phases (closure phases). The 2011 data set benefited Note-Here we list the UT date, the average modified Julian date of the night of observation, the baselines used in their corresponding configuration, the number of useful squared visibility points obtained for the night, the number of useful closure phase points obtained for the night, the rotation phase for the primary star in λ And, and the calibrator stars that were used for each corresponding night. The rotation phase is derived by using the first observation in 2010 as the zero point.
from MIRC having been upgraded earlier that year, allowing for simultaneous use of all six telescopes. These upgrades provided data sets to acquire up to 15 visibilities, 10 independent triple amplitudes, and 10 independent closure phases for each spectral channel. Parks et al. (2021) detail the reduction steps and error corrections but we will briefly note some of their steps here. The data were reduced using the official IDL pipeline for reducing MIRC data (Monnier et al. 2007). Each block of raw fringe data contained coadded frames, and were corrected for any instrumental effects by background subtraction in order remove instrumental noise and foreground normalization to correct for any pixel-to-pixel variation. Raw square visibilities, closure phases, and triple amplitudes are output through the use of Fourier transforms and are photometrically calibrated. The data were corrected for the atmospheric coherence time and optical changes in the beam path with the use of calibrator stars that were taken either immediately before or after the target λ And.

Data Reduction
In the 2010 data, one of the calibrators 37 And (HD 5448) was found to be a binary by Che et al. (2012) and had its orbit fully characterized by Roettenbacher et al. (2016b). Parks et al. (2021) formed a comparison of using 37 And as either a single star calibrator or as a binary calibrator. They found that these comparisons only incurred an error of 1.24% for the square visibilities, which is well below the multiplicative error correction, and a closure phase standard deviation of 1.14 • . We ex- The differences in angular sizes for 7 And and σ Cyg used between the two years are small and within their 1σ errors. a This is the semi-major axis angular separation of the binary calculated by Roettenbacher et al. (2016b). ecute a separate reduction and calibration for the 2010 data set using the official MIRC reduction pipeline in order to correct for the 37 And binary calibrator. We use the more recent calibrator diameter estimates, whose values differ from Parks et al. (2021), for this new reduction and calibration. The data uncertainties also go through a post-calibration process to account for known systematic errors of the MIRC instrument.
For the 2010 data we kept the same systematic errors as Parks et al. (2021). These errors are different compared to the 2011 data set as the quality of the 2010 data are taken with a four telescope configuration and are of lower quality while the higher quality 2011 data are taken with a six telescope configuration. A 15% multiplicative error correction was used in association with the transfer function, a 2 × 10 −4 additive error correction was used in association with bias at low amplitudes for the square visibilities, and a 20% multiplicative error correction and a 1 × 10 −5 additive error correction was used for the triple amplitudes. The same 1 • error floor was used for the closure phases as was used in Zhao et al. (2011). We present the square visibilities and closure phases for the 2010 data set in Figure 1.
We use the same calibrator diameter estimates listed in Parks et al. (2021) since the 2011 data set has been reduced and calibrated. Even though different angular sizes were used for the calibration of the 2010 and 2011 data set for 7 And and σ Cyg, the differences between the two angular sizes reported in Table 2 are small and within their respective 1σ errors. Systematic errors were taken into account during calibration similar to that of Monnier et al. (2012). A 10% multiplicative error cor-rection was used in association with the transfer function for the 2011 data and a 2 × 10 −4 additive error correction was used for the square visibilities. A 15% multiplicative error correction was used and a 1 × 10 −5 additive error correction was used for all the triple amplitude data. Lastly, the same 1 • error floor was used for the closure phases just as it was presented in Zhao et al. (2011). We present the square visibilities and closure phases for all of the 2011 data set in Figure 2.

MODELING λ AND WITH SIMTOI
The SImulation and Modeling Tool for Optical Inteferometry (SIMTOI) is an interferometric modeling code 1 (Kloppenborg & Baron 2012a,b;Kloppenborg et al. 2015) that uses a Graphical Processor Unit (GPU) to represent stars and their environments in a threedimensional framework. In SIMTOI, the stellar intensity maps are two-dimensional textures applied on top of orbiting/rotating three-dimensional stars. Once the scene is rendered, the GPU also powers the fast computation of interferometric observables. SIMTOI offers a large choice of global and local optimizers to solve Maximum A Posteriori (MAP) or model selection problems. Our first goal in using SIMTOI was to derive initial guesses for λ And's stellar parameters (such as its rotation axis), since our imaging code would be too slow to wade through the entire parameter space. Our second goal was to assess the potential number of spots present on the star via model selection. Both tasks were solved using the MultiNest optimizer (Feroz & Hobson 2008;Feroz et al. 2009Feroz et al. , 2019, which implements the Importance Nested Sampling algorithm.

Modeling λ And
We devised models of λ And with different number of circular spots, from three to six. Six parameters were used to model the star itself: rotation period, rotation axis (inclination and position angles), temperature, angular diameter, and coefficient of the power limb-darkening law (Hestroffer 1997). The stellar parameters were given uniform prior distributions within a wide range of values, based on the stellar parameters listed in Parks et al. (2021) as a starting point (e.g. ±20 • for angular parameters). Four parameters were used per spot: longitude, latitude, diameter, and flux. These spot parameters were also given uniform distribution. In particular, their location was not constrained.
For each data set -2010 or 2011 -SIMTOI renders an image per epoch (day). The rendering resolution was set to a 64 × 64 image with a 0.05 mas/pixel resolution. MultiNest was run for each model and converged after a few hours, providing MAP parameter values, as well as the marginal likelihood values (the so-called log Z).

Modeling results
We report the χ 2 and log Z values for each spot model in Table 3. We also provide the approximate nominal values for the physical parameters. MultiNest does provide error bars, but since they do not account for sys- tematic errors, they are vastly underestimated. While one could bootstrap the data before MultiNest runs, this would be too computationally intensive and yet still imprecise due to our approximate modeling of spots. Our model spots are circular, which may be an unrealistic assumption, but is sufficient to identify the main potential location of intensity peaks on the surface. The log Z values are maximal for the five spot model for the 2010 data, and the four spot model for the 2011 data. The corresponding reduced χ 2 values are low for the 2011 data and much higher for 2010. Setting aside the possible differences in error calibration between 2010 and 2011, this would indicate that the 2010 surface map is much more complex than the 2011 one (which we did confirm during imaging).
We ultimately choose the 4 spot model for the 2011 data as the best representative model that produces the most accurate parameterization of λ And. The estimated 54.2 day rotation period of the primary from our model using the 2011 data set is consistent with other works. Henry et al. (1995) reports a rotation period 54.07 days from their photometric analysis while Parks et al. (2021) reports a 54.02 ± 0.88 day rotation period from their own photometric analysis and an average of a 56.9 ± 8.8 day rotation period from their interferometric analysis. While the 2010 data set had a larger rotation phase coverage than the 2011 data set, the rotation period based on the 2011 data are overall more reliable based on MultiNest results and the fitting of the model to the data. This is most likely due to the larger Note-Higher ln Z value is better, lower χ 2 ν is better. No error bars are calculated since the models from SIMTOI using Multinest does not currently generate reliable error bars. We rely on the imaging results for more precise measurements and calculation of errors.
amount of (u, v) coverage, number (u, v) points, triple amplitudes, and closure phase points in the 2011 data set compared to the 2010 data set. This calculated period from the four spot model using the 2011 data is consistent with previous works.  Baron & Martinez in prep) which models the stellar surface temperatures of single stars or binary systems as two-dimensional arrays on top of a stellar geometry. The stellar geometry itself is defined either by analytic formulas (ellipsoids, fast rotators) or by solving Roche equations. In imaging and model-fitting problems, ROTIR makes use of the optimization packages OptimPack (Thiebaut 2002) and NLopt (Johnson 2008) to maximize the posterior probability of the model.

Geometry setup
ROTIR relies on the package OITOOLS 3 ( Baron et al. 2019) to read in our data, split up or combine our data temporally, and plot any images featured in this work.
Once the interferometric data are read, we define the stellar parameters and orientation of our object. Our code requires several parameters: the angular size at 2 https://github.com/fabienbaron/ROTIR.jl 3 https://github.com/fabienbaron/OITOOLS.jl the pole in milliarcseconds, the surface temperature, the fraction of the critical angular velocity if the star is rapidly rotating, the limb-darkening law and its corresponding coefficient(s), the exponent needed if there is any gravity darkening (von Zeipel 1924), the difference in angular velocity between the equator and the pole, the inclination, position angle, and rotation period of the star. Our code allows the user to choose between three different limb-darkening laws: a quadratic law, logarithmic law, or Hestroffer law (commonly known as the power law; Hestroffer 1997).
Our geometrical setup starts with selecting a tessellation scheme. Two schemes have been implemented so far: the HEALPix tessellation (Górski et al. 2005) and the latitude/longitudinal scheme. HEALPix presents the advantage of equal area tessels, provided the star does not depart too much from a spherical shape. The latitude/longitudinal scheme allows for simulating differential rotation, but requires more tessels to represent the surface. As part of this work we tested both tessellation schemes, which result in qualitatively identical maps. Most results presented in this paper were obtained with the latitude/longitude scheme. The number of pixels per angular diameter was chosen based on the estimated angular diameter size divided by the imaging resolution limit. Therefore, the minimum total number of pixels required across the surface of a star would simply be the number of latitude pixels times the number longitude pixels.
For the latitude/longitude scheme, the number of latitude pixels is based on the number of pixels per angular diameter since the latitude range spans from −90 • to 90 • and number of longitude pixels is twice the number of pixels per angular diameter since the longitude ranges from 0 • to 360 • . We number the vertices of the polygon by 1, 2, 3, 4 in a counterclockwise direction when viewed along the direction of the normal. A fifth element is also included for each pixel and defined to be at the center of each pixel.
Once the user has chosen a tessellation scheme and calculated the number of pixels required for imaging, the user then has the choice of choosing between three different geometries: a scaled unit-sphere, an oblate spheroid, or a Roche object 4 . The order in which the pixels are mapped out on the surface of the star are counterclockwise when viewed along normal of the positive z direction on the (x, y, z) plane.
The surface area A n is calculated for all n pixels in order to determine the amount of relative flux coming from the star with the following where v is the vector of (x, y) projected positions of the n th pixel in a 2-dimensional (x, y) plane at the j th corner with m number of corners in the polygon of choice, · is the scalar product, and ∧ the vector cross product operator. The m + 1 corner here points back to the first corner of the pixel.
Once the surface area of each pixel is calculated with the desired limb-darkening law, the Fourier transform S is done on every pixel for a 3-dimensional object (Lee & Mittra 1983;Chu & Huang 1989;McInturff & Simon 1991) in order to compare the frequencies of our data on the (u, v) plane by using the following equation where k is a vector containing each u and v frequency on the (u, v) Fourier plane. We use the flux to visibility matrix S to compute the model visibilities using: where V is the model complex visibility vector, T is the temperature map vector, L is the limb-darkening map, • is the Hadamard (element by element) vector product, and the division is the Hadamard division.

Differential rotation option
The user can select whether or not to turn on the option to simulate differential rotation. The equation for differential rotation  used in our code is in the form where Ψ is the latitude, Ω(Ψ) is the rotation rate at a specific latitude, Ω eq is the rotation rate at the equator, and ∆Ω is the difference in angular velocity between the equator and the pole. This difference between angular velocity in the equator and the pole is related to the differential rotation coefficient, k, or the surface shear parameter, α, commonly found in the literature (e.g., Henry et al. 1995;Davenport et al. 2015;Kővári et al. 2015) and is defined through the following equation or in terms of the polar and equatorial rotational periods as k = P pole − P eq P pole .

First use of regularization
Fitting a model to the data with no prior constraints will produce unrealistic images due to overfitting. The Maximum A Posteriori method balances the likelihood term with our prior expectations of what the temperature map should look like. The optimal temperature map x opt is then found as: where χ 2 (x) is the chi square fit of the data to the model, R(x) is the regularization, and µ is the hyperparameter setting the relative weight of the regularizer versus the likelihood. We implemented three different regularizations for use in ROTIR: positivity, l 2 norm, and total variation. Positivity enforces a non-negative temperature map. The l 2 norm takes the square root of the sum of square values for each pixel and penalizes pixel values straying too far from the average value. Our third regularizer is total variation which computes the spatial gradient of the model image and penalizes large temperature fluctuations between neighboring pixels such that it shows smoother transitions on a local scale.

APPLYING ROTIR TO λ AND
For λ And, we use positivity and total variation as the two regularizers necessary to determine the best image. Using the l-curve method (Renard et al. 2011), we choose a weight of µ = 0.01 that has a small amount of regularization before entering into a regularization dominated regime. We show examples of strong and weak regularization in Figure 3 for data of λ And taken on 2011-Sep-14 to prove why we need a good balance between regularization and pure model fitting when finding an optimum image.
In order to determine of the number of tessels required on the surface of the star, we use the parameters we obtained from modeling λ And using SIMTOI. Knowing that the CHARA angular resolution limit is θ ≈ 0.60 mas at H-band (λ = 1.61µm), we estimate that we need 40 pixels across the visible equator to meet Nyquist sampling (imaging resolution limit is θ ≈ 0.30 mas in H-band). Therefore, we use 80 pixels around each latitude, including pixels behind the star, and 40 pixels across each longitude for a total of 3200 pixels on the surface of the star. Our sampling of pixels across the resulting images are solely based on the number of pixels on the surface on the star and not the overall field as the field size can be arbitrarily chosen based on the plotting axes.

A first look at imaging
In order to find the best geometrical setup for primary star in λ And, we test both a spherical star and a Roche lobe shape to see if there is any signs of major Roche lobe overflow. While Donati et al. (1995) and Parks et al. (2021) both suggest that there is no Roche lobe overflow, we decide to investigate this for λ And since slight oblateness was found in another RS CVn variable, ζ Andromedae (Roettenbacher et al. 2016b).
We start with the parameters from SIMTOI to create our spherical star and, with the addition of the longitude of the ascending node, argument of periapsis, and eccentricity found in Walker (1944), create our Roche lobe geometry. Donati et al. (1995) states that λ And is coplanar, therefore we use the inclination rotation axis of the primary star as the inclination of the orbit. We use the same hyperparameter and apply a uniform temperature map across the whole star as an initial condition for both geometries. Using a Julia package called OptimPack 5 that solves for an optimum temperature map through a quasi-Newtonian method (Thiebaut 2002), we obtain for the best temperature map given our all our data in a 5 https://github.com/emmt/OptimPackNextGen.jl given year. This algorithm compares the Fourier transforms from Section 4.1 to the 2011 data to solve for the best temperature map.
The resulting criterion for the Roche lobe geometry is higher (χ 2 (x)+µR(x) = 6288) when directly comparing it to a spherical geometry (criterion = 4489). We also find that the pole-to-equator ratio at the L1 Lagrangian point for the primary is 0.9967. With these two calculations, we determine that a spherical geometrical shape for the primary of λ And is an acceptable approximation for the true shape of the star.
Once we have determined that the spherical geometrical setup is the most optimal for λ And and choose the most optimal regularization weight, we are now set for calculating the best fit for the temperature map. We present the resulting Mollweide maps of λ And for both epochs in Figure 4. These maps reflect no time variability and assume that λ And is undergoing solid-body rotation. A better representation of the temperature maps are shown in Figure 5 for each given night in 2010 and 2011.
A first look at the temperature maps between the 2010 and 2011 epochs shows a few interesting characteristics about λ And's surface. Comparing the two temperature maps show notable similarities for two spots in the northern hemisphere between the two epochs (i.e., the spot around 20 • latitude and ∼ 100 • longitude, and the spot around 0 • latitude and 170 • in both epochs). There are two other notable spots that either disappear or appear from one epoch to the next. The spot in the 2010 epoch around 30 • latitude and 150 • longitude seems to has disappeared within the 2011 epoch. A spot seems to be forming within the 2011 map in the southern hemisphere around −40 • latitude and 50 • longitude with hints of its emergence with similar place in the 2010 epoch. We note that the spot in the 2010 epoch around 15 • latitude and −90 • does not appear in the 2011 epoch. This is most likely due to missing rotational phase coverage in the 2011 data set.

Refinement of physical parameters
After finding the best model from SIMTOI, we use the parameters from the 4-spot model based on the 2011 data and use the bootstrap method. We apply the bootstrap method in order to find the final parameters and errors for the primary component of λ And. We use 50 bootstrap iterations to solve for only four parameters: angular radius, the limb-darkening coefficient, inclination, and position angle. We choose to leave the rotation period of the primary fixed throughout this bootstrap because there is a degeneracy towards lower rotation periods. We believe that this is due to the fitting algo-  rithm in ROTIR choosing the difference between the first and last observing date (in a given epoch) for the period instead of the true rotational period. Our bootstrap is dependent on the NLopt package (Johnson 2008) and Nelder-Mead Simplex method (Nelder & Mead 1965;Box 1965;Richardson & Kuester 1973) within NLopt for obtaining our final parameters with their corresponding errors. We restrict lower and upper bounds within NLopt for these four parameters as follows: their associated errors are calculated though their standard deviation. Our 50 bootstraps do not show a Gaussian distribution, but we are prevented from running a large number of bootstraps due to computation time. We indicate that parameters do not deviate too largely from their mean values. It is likely that doing more bootstraps will slightly increase the error bars but not in a significant manner. We show the results of our bootstrap values in Figure 6.

Images of λ And
Temperature maps are not indicative of what is actually represented from observations. In order to present  an image, we include to use a power law for limbdarkening (Hestroffer 1997) and multiply it by the cells of the temperature maps that are visible to the observer. We use the limb-darkening coefficient from our bootstrap to present the images in Figure 7 and present the physical parameters for the primary star in λ And using the parameters from our bootstrap in Table 4. Here we compare images made independently from ROTIR to another image reconstruction code called SURFace imagING (SURFING) in Figure 8. SURFING is a Monte Carlo based imaging code written in IDL specifically written for imaging spheroids (see Roettenbacher et al. 2016b). Overall, there is a good agreement between the two imaging methods. Since we are only focusing on the imaging comparison aspect for these two codes, we see that the spot locations and contrast between the two are very similar, with a few minor differences, as shown in Figures 5 and 8.

Comparison to Parks et al.
The results of this work largely agree to those of Parks et al. (2021) with exception of the inclination of λ And being the only disagreement. Parks et al. (2021) used a combination of a genetic algorithm (Charbonneau 1995) and the Nelder-Mead Simplex method (Nelder & Mead 1965;Box 1965;Richardson & Kuester 1973), in order to make individual models for each night of data. Each surface model calculates an angular diame-   Figure 5) were made by using two different consecutive nights and merging the data as one night. We find that this does not largely affect the results of the imaging since the rotation made from two consecutive nights only span ∼2% of the rotation period.
ter, limb-darkening coefficient based on the power law, a starspot covering factor, starspot latitude, starspot longitude, and starspot intensity ratio for λ And. Once all the models were made, Parks et al. (2021) traced each starspot on the surface for each epoch. Ellipse fits to starspot positions were calculated, and an average computed position angle and inclination angle were made from these ellipse fits for each year. Parks et al. (2021) reported that the inclination of primary from their 2010 and 2011 data is 75 ± 5.0 • and 66.4 ± 8.0 • , for each respective year, giving an overall average of 70.35 ± 6.7 • while we report an inclination of 85.63 ± 2.32 • . We believe that our calculations from this work are accurate for several reasons. The ini-tial SIMTOI calculations were done with a global search with no restrictions in parameter space, including inclination. The resulting parameters obtain from SIMTOI were then used in ROTIR with a sufficient range that included the inclination value from Parks et al. (2021). If the value for our inclination were incorrect and actually leaned towards this previous value, the resulting bootstrap method would have reflected it by converging on the lower bounds of our parameter space using our bootstraps. In addition, the work by Parks et al. (2021) relied on independent models for each night and tied them together to form an analysis while we use the all the data of each epoch collectively to form one image.  b Based on the physical radius from this work and the log g from Drake et al. (2011).
c Based on the physical radius from this work and the effective temperature from Drake et al. (2011).
d Since Parks et al. (2021) had multiple values reported for the same parameter, we show the averages of the respective parameter here.
Note-The observed parameters were optimized through a bootstrap approach with the exception of the rotation period, which was fixed. We take our fixed rotation period parameter directly from the best model in SIMTOI.
7. BEYOND SOLID ROTATION IMAGING 7.1. Simulating differential rotation In our Figures 4 -8 using SIMTOI/ROTIR, SURFING, and in Parks et al. (2021), all imaging has been performed assuming that the star is rotating as a solid body however, we attempt to estimate differential rotation through our data. Henry et al. (1995) studied photometry of λ And over 14 years and found evidence of shear across the surface. In order to see if we are able to detect any differential rotation with our interferometric data, we simulate starspots on a star with a low differential rotation coefficient and a low temperature gradient on the surface, and have the spot move across a few days with the same period as λ And. Then we do a cross-correlation for each latitude band on the star and see if there is any deviation from zero.
Our simulations show two different scenarios. The first simulation presents a highly unrealistic starspot that are two pixels wide in longitude and spanning throughout all latitude from pole to pole. Our second simulation shows two circular starspots that are 5 pixels in radius at +45 • and −45 • latitude (in respect to the equator) and at 135 • longitude. We presents our simulations of a simple star with similar parameters as λ And using differential rotation coefficient from Henry et al. (1995) of k = 0.04, which corresponds to differential angular velocity (∆Ω) of 0.26, in Figure 9.

Testing differential rotation on λ And
We apply the same cross-correlation method for the 2011 data set and calculate the deviations. We find that we are unable to detect any differential rotation with our data due to three reasons. First, our data does not span an entire rotation, therefore we are not able to compare the same spots from the previous rotation. Second, λ And is a very slow rotator so we do not have enough resolution to detect any small amounts of differential rotation, if differential rotation truly exists on λ And. In fact, the large scale magnetic spots on λ And may not be able to be used to measure any real surface differential rotation based on its dynamo. Korhonen & Elstner (2011) states that surface differential rotation can only be recovered by observing the spot motion of small spots, unlike λ And's large scale magnetic spot structure. Third, the amount of square visibilities and closure phases for each observation are sparse for most observations. Since the goal is to detect any shear as evidence for differential rotation, we reconstruct an individual temperature map for each observation date from Figure 9. We show simulations of differential rotation by doing a correlation using the unrealistic starspot among a longitudinal band (left) and two starspots (middle) of a fake star with the same parameters of λ Andromedae within the 2011 epoch (with the exception of the temperature map). The differential rotation coefficient we use here is ∆Ω of 0.26 from Henry et al. (1995).
The plot (right) shows the number of pixels that have shifted in respect to the longitude after subtracting off the total shift of a spot for a given latitude. The pink line at this coefficient represents the unrealistic starspot change in pixels while the yellow line shows the two starspots change in pixels as a function of the longitude. We choose to compare the first and last observations within the 2011 epoch to show the maximum amount of correlation.
the 2011 epoch but initialize with the temperature map obtained from Figure 4. We show our results in Figure  10.

Updated Orbital Parameters and Secondary Parameters
Using the updated parameters from the primary star in λ And in this work, the mass ratio from Donati et al. (1995) and Kepler's Third Law, we are now able to calculate the mass of the secondary and the semi-major axis of the binary system. We calculate that the mass for the companion is 0.15 +0.09 −0.05 M and with corresponding semi-major axis is 6.12 mas for the system.

The Search for the Secondary
We begin our search for the companion by obtaining an estimate on the luminosity ratio and angular size of the secondary to narrow down our search. For the luminosity ratio, we used a mass-luminosity relation for each corresponding star in our system (L 2 /L 1 = 0.23(M 2.3 2 /M 4 1 )) and calculated to be approximately L 2 /L 1 = 0.00121. If we assume that the H-band flux ratio is the same as the luminosity ratio of the two stars and if we use the H-band magnitude of the primary 1.40 mag (Ducati 2002), this would correspond to an estimated H-band magnitude of 8.7 mag for the secondary. This is slightly beyond MIRC's magnitude limit and not likely to be detected, however we still investigate the possibility of detection. In order to calculate the estimated angular size of the secondary, we first calculate the physical size by using the mass-radius relation (R = 0.0753 + 0.7009M + 0.2356M 2 ) developed by Maldonado et al. (2015) for low-mass stars. Given that the calculated physical radius is 0.19 R , we find that the estimated angular radius is approximately 0.03 mas. Now that we have an estimation of the angular size and flux ratio, we perform a grid search in right ascension and declination over a 10 mas distance from the primary star for every night in the 2011 epoch. This approach is similar the methods used in Baron et al. (2012) and CANDID (Gallenne et al. 2015) with the difference that the primary is using the model visibilities obtained during image reconstruction. We model binary visibilities and vary both the brightness ratio and the angular radius for the secondary using NLopt for each section of the grid. We restrict the parameter space for the angular radius to [0.0, 1.0] mas while restricting the flux ratio (secondary/primary flux) for the system from [0.0, 0.2].
While we do find that the average flux ratio using the 2011 data set of 0.00213 ± 0.00116 is within the theoretical estimated value, we find two major reasons for believing that we were not able to find the secondary companion. First, the average angular radius found by using the 2011 data is 0.602±0.356 mas, largely inconsistent with our estimation using mass-radius relation for low-mass stars. Our errors for both the flux ratio of the system and the angular size of the secondary were calculated by taking the standard deviation of every night's grid search result from the 2011 epoch. The values of angular radius for an individual night were also seen to hit a boundary condition (either 0 mas or 1 mas), thus assessing that the calculated values are incorrect. Second, the best fit right ascension and declination positions for each night in the 2011 data set were positioned in a random assortment on the grid space with no clear indication of a circular or elliptical orbit.
Another reason that we may not be able to find the secondary for λ And could be due to lack of (u, v) coverage for each individual night in the 2011 epoch data set. For this reason, we proceed to not use the 2010 data set to find the secondary as those observations were taken with two different sets of 4T observations in a given night and as a result do not provide better (u, v) coverage compared to the 2011 data set.

CONCLUSION AND FUTURE PROSPECTS
In this paper, we do interferometric modeling and imaging on λ And for the 2010 and 2011 epochs. First, we use SIMTOI in order to find which model is most probable for finding the best parameters. Then we use the parameters from SIMTOI and use them for imaging in ROTIR. Using the parameters from the best SIMTOI model as a starting point, we apply the bootstrap method to get the final physical parameters for λ And. We find that our images from ROTIR fairly agree with the images produced to the other image reconstruction code, SURFING, and our physical parameters are also fairly consistent of previous works with the exception of the inclination.
Images from both codes show that the spots on λ And from both epochs seem to favor certain latitudes and are mostly concentrated in the northern hemisphere. For both the 2010 and 2011 epochs, we find that most of the spots are centered around +20 • latitude. These spot concentrations to a certain latitude are consistent with the interferometric images shown in Roettenbacher et al. (2016b) of ζ Andromedae, another RS CVn variable. The absence of symmetrical spots on active latitudes as observed on the Sun is evidence that λ And may not have a solar-like dynamo.
Finally, once we produce static images of the primary star in the system, we test to see if we find any evidence for differential rotation and detect the secondary companion. We start with a simulation of differential rotation and compare those results to the 2011 interferometric data set. Our results remain inconclusive as we cannot detect any sheer within the 2011 data set largely due to λ And being a slow rotator. In our search for the companion, we do a grid search by fitting various models for the companion (i.e., varying the angular radius of the secondary and flux ratio of the system). While the flux ratio was consistent with the approximated value, the angular radius was largely inconsistent with our estimated calculation therefore concluding that we were unable to detect the secondary.
Our ROTIR code is not just limited to interferometric imaging but is also capable of light-curve inversion. Our future work will plan on using the multi-band photometry in Parks et al. (2021) and compare those resulting images with the interferometric images from this work. Our plans also include using the photometric data as a bridge for the 2010 and 2011 interferometric epochs in order to detail how λ And is evolving over the course of a year. We are currently implementing additional numerical techniques to ROTIR (Abbott et al. in prep) in order to improve light-curve inversion quality with the use of Alternating Direction Method of Multipliers (Chan et al. 2011). Finally, we have future plans to implement Doppler imaging and Zeeman-Doppler imaging into ROTIR.