Tackling the Unique Challenges of Low-frequency Solar Polarimetry with the Square Kilometre Array Low Precursor: The Algorithm

Coronal magnetic fields are well known to be one of the crucial parameters defining coronal physics and space weather. However, measuring the global coronal magnetic fields remains challenging. The polarization properties of coronal radio emissions are sensitive to coronal magnetic fields. While they can prove to be useful probes of coronal and heliospheric magnetic fields, their usage has been limited by technical and algorithmic challenges. We present a robust algorithm for precise polarization calibration and imaging of low-radio frequency solar observations and demonstrate it on data from the Murchison Widefield Array, a Square Kilometer Array (SKA) precursor. This algorithm is based on the Measurement Equation framework, which forms the basis of all modern radio interferometric calibration and imaging. It delivers high dynamic range and fidelity full Stokes solar radio images with instrumental polarization leakages $<1\%$, on par with general astronomical radio imaging, and represents the state-of-the-art. Opening up this rewarding, yet unexplored, phase space will enable multiple novel science investigations and offer considerable discovery potential. Examples include detection of low-level of circularly polarization from thermal coronal emission to estimate large-scale quiescent coronal fields; polarization of faint gyrosynchrotron emissions from coronal mass ejections for robust estimation of plasma parameters; and detection of the first-ever linear polarization at these frequencies. This method has been developed with the SKA in mind and will enable a new era of high fidelity spectro-polarimetric snapshot solar imaging at low-radio frequencies.


INTRODUCTION
The solar corona is the outermost layer of the Sun and is made up of hot magnetized plasma. The magnetic field permeating this plasma couples the solar atmosphere to the solar interior. Solar magnetic fields are responsible for the bulk of the observed coronal phenomenon, spanning a range of time scales from a solar cycle to flares lasting milliseconds and in terms of the energetics from massive coronal mass ejections (CMEs) to nanoflares. Coronal magnetic fields also play a large role in shaping the coronal structures and even the solar wind. Hence, to understand coronal dynamics it is essential to measure and understand the ever-evolving coronal magnetic fields.
The sun and the corona is routinely observed using multiple space-based and ground-based instruments, spanning from radio, visible, extreme ultraviolet (EUV), to X-ray wavelengths. Magnetic fields, however, are rather hard to measure (Wiegelmann et al. 2014). Routine measurements of global solar magnetic field are available only at the photosphere (Lagg et al. 2017). For lack of better options, large-scale quiescent coronal magnetic fields can only be estimated using extrapolations of photospheric magnetic fields (e.g. Wiegelmann 2008;Jiang et al. 2014). During CME eruptions, EUV (Gopalswamy et al. 2011, etc.) and white light (Schmidt et al. 2016, etc.) have also been used to measure the coronal magnetic fields. Very recently, the global coronal magnetic field has been measured successfully using nearinfrared observations in the range 1.05-1.35 R , where R is the solar radius (Yang et al. 2020). Measuring coronal magnetic fields is also a key science objective for the Daniel K. Inouye Solar Telescope (DKIST; Rast et al. 2021), though these measurements are likely to remain limited to < 1.5 R .
Under favorable circumstances, radio observations have also been used to estimate coronal magnetic fields associated with active regions and/or CMEs (Alissandrakis & Gary 2021). The coronal optical depth at these heights (> 1.3 R ) becomes too low for visible and EUV bands. Although the radio observables, in principle, are sensitive to coronal magnetic fields, it has been technically too challenging to extract this information on a regular basis. Most of the radio studies have focused on the active emissions, and the global quiescent coronal magnetic fields at high coronal heights have remained beyond reach.
The polarization properties of solar radio emission, in addition to being a direct probe of the coronal magnetic fields, can also provide strong constraints on the emission mechanisms. Despite its well-appreciated importance, low-frequency polarimetric observations of the Sun are one of the least explored area of solar physics. The radio Sun is a complicated source. It has structures spanning a large range of angular scales. The spectral, temporal, and morphological characteristics of the radio emissions are also very dynamic. The brightness temperature (T B ) of the low-frequency solar emissions can vary from 10 3 K for gyrosynchrotron emission from CME plasma (e.g. Bastian et al. 2001;Mondal et al. 2020b) to 10 13 K for bright type III radio bursts (e.g. McLean & Labrum 1985;Reid & Ratcliffe 2014) over a background quiescent T B of ∼ 10 6 K. Depending upon the emission mechanism at play, the polarization fraction can vary from ≤ 1% to ∼ 100% (e.g. McLean & Labrum 1985;Nindos 2020).
To date, most of the polarimetric studies of the Sun at low frequencies are based on nonimaging dynamic spectra measuring the circular polarization (e.g. Reid & Ratcliffe 2014;Kaneda et al. 2017). These observations cannot provide any information about the source structure or location. Some innovative instruments use simultaneous Stokes I imaging and Stokes V dynamic spectra to help in localization of the source of active emission (Sasikumar Raja 2014). These studies implic-itly assume the locations of the peaks in the total intensity and circularly polarized emission to be the same. This assumption usually holds when there is a single dominant source of emission. This approach is not useful when multiple sources of active emission are simultaneously present on the Sun (Mohan & Oberoi 2017) or for weaker and/ or extended emission like gyrosynchrotron emission from CME plasma (e.g. Bastian et al. 2001;Mondal et al. 2020b) and the free-free emission from the quiet Sun (Sastry 2009).
The variation in the solar emission over small temporal and spectral scales imposes a requirement for snapshot spectroscopic imaging. The need to be able to see features varying vastly in T B highlights the need for a high imaging dynamic range. Only recently it has become possible to meet these exacting requirements for solar radio imaging. Several new generation low-frequency radio interferometers have now become available -LOw Frequency ARray (LOFAR; van Haarlem et al. 2013) operating at 10 -80 MHz and 120 -240 MHz; Long Wavelength Array (LWA; Kassim et al. 2010) at 10 -88 MHz; and the Murchison Widefield Array (MWA; Tingay et al. 2013;Wayth et al. 2018) at 80 -300 MHz. Of these, the MWA has a large number of antenna elements (a total of 128 antenna elements) spread over a small footprint in a centrally condensed configuration and is especially well suited for high-dynamic-range high-fidelity snapshot spectroscopic imaging. Even with the availability of data from a capable instrument like the MWA, there remains another problem to overcome -that of being able to produce the large number of high-quality images necessary to enable the science in a reliable unsupervised manner. This need was met by a robust Stokes I calibration and imaging pipeline developed by Mondal et al. (2019) designed to build on the strengths of the MWA architecture. This pipeline is christened Automated Imaging Routine for Compact Arrays of the Radio Sun (AIRCARS) and the highfidelity spectroscopic snapshot images it delivers now represent the state of the art and have already led to new discoveries such as quasiperiodic pulsations of solar radio bursts (Mohan et al. 2019b,a;Mondal & Oberoi 2021), weak nonthermal emissions from the quiet Sun (Mondal et al. 2020a;Mondal 2021;Mohan 2021a,b), and weak gyrosynchrotron emission from CMEs (Mondal et al. 2020b). Having established the ability of MWA solar observations to deliver high-fidelity Stokes I images, it is now the right time to use that experience for high-fidelity polarization imaging.
Full-Stokes calibration is significantly more challenging than working with Stokes I alone. These challenges are even greater for the case of low radio frequency so-lar imaging. On the one hand, solar emission can have a very large range of intrinsic polarizations, which can also vary rapidly across time and frequency. The MWA has a large field of view (FoV). Based on the FWHM of the primary beam, at 150 MHz the FoV of the MWA is ∼610 degree 2 , which reduces to ∼375 degree 2 by 200 MHz (Tingay et al. 2013). The wide FoV aperture arrays tend to have large instrumental polarization imposing a strong requirement for precise calibration. In fact, some of the assumptions made for routine polarimetric calibration at higher frequencies for small FoV instruments no longer hold in this regime (Lenc et al. 2017). We present a general algorithm for polarimetric calibration suitable for our application. We have implemented this algorithm and demonstrate its efficacy on the MWA solar data. The algorithm will be well suited for solar imaging with the future SKA-Low and other interferometers with centrally condensed array configurations.
The paper is organized as follows. We first briefly discuss some basics of polarization calibration in Section 2 to build up the base for our calibration algorithm. Section 3 describes the challenges of the polarization calibration of the Sun at low frequencies and the limitations of the conventional methods of polarization calibration. We then describe the new algorithm in Section 4. Section 5 presents the final results with a discussion, and Section 6 provides the conclusion.

POLARIZATION CALIBRATION OF A RADIO INTERFEROMETER
A radio interferometer is made up of a number of radio antennas or elements. These antennas measure the voltages corresponding to the two orthogonal polarizations of the electric field, E, incident on the antenna. E could be measured in either of the linear or circular bases -(E X , E Y ) or (E R , E L ) respectively. Incident E induces a voltage in the antenna, and the primary observable of an interferometer is the cross correlation between the components of the induced voltages for every antenna pair, referred to as visibilities. To capture the complete information about the state of polarization of E, for any given antenna pair described by indices i and j, an interferometer needs to measure X This complete set of visibilities is often referred to as full-polar visibilities.
The measured visibilities include the corruption due to atmospheric propagation effects and instrumental effects. In order to arrive at the true visibilities corre-1 † represents conjugate transpose.
sponding to the astronomical sources, these corruptions need to be removed. This process is known as calibration.  proposed a general framework for polarimetric calibration based on the measurement equation. Briefly, the measured complex voltage vector, V, per antenna can be expressed in terms of E, and the antenna-based Jones matrix, J, (Jones 1941) as: where J is a 2 × 2 matrix representing the instrumental and atmospheric effects. The measured correlation products for two antennas represented by indices i and j can be written as a 2 × 2 matrix, also known as the Visibility matrix, as follows (Hamaker, J. P. 2000;Smirnov 2011): where V ij is the true source visibility matrix; V ij is the observed visibility matrix; and V I , V Q , V U and V V are the Stokes visibilities of the incident radiation. Here we have followed the IAU/IEEE definition of the Stokes parameters (IAU 1973;. The four Stokes parameters, I, Q, U, V were originally defined by Stokes (1851). Stokes I represents the total intensity; Stokes Q and U represents the linear polarization, and circular polarization is denoted by Stokes V. All V ij s have independent additive noise, N ij , associated with them, and the Eq. 2 can be written as: Equation 3 is referred to in literature as the measurement equation of a radio interferometer Hamaker, J. P. 2000;Smirnov 2011). The objective of polarization calibration is to estimate J i s for all antennas and obtain V ij from the V ij . This requires correcting for four different aspects. These aspects and their impacts are enumerated below.
1. Time-variable instrumental gain: While they are independent in origin, in practice, it is not feasible to disentangle the atmospheric propagation effects from time-variable instrument gains. So these effects are clubbed with instrumental gains in this formalism. The impact of these gain variations is to make the interferometer incoherent.
2. Frequency-dependent instrumental bandpass: This can modify the true spectral signature of the source and introduce incoherence across the frequency axis.
3. Polconversion: This can lead to a leakage from Stokes I to other Stokes parameters and thus modifies the observed magnitude of the observed polarization vector, p = (Q, U, V).
4. Polrotation: This is the mixing between Stokes Q, U, and V and leads to a rotation of p.
The time and frequency dependence of the instrumental gains arise due to the nature of the signal chain and the atmospheric propagation effects. These are routinely corrected for in standard interferometric calibration. Ideally a set of orthogonal receptors are expected to receive only the matched orthogonal component of the incident E. In practice, reasons ranging from proximity to other receptors, manufacturing tolerances to cross talk between closely placed cables, imply that a given orthogonal receptor also picks up some amount of signal of the other component of E. This mixing of the orthogonal components of E gives rise to polconversion. Things like misalignment of the dipoles with respect to the sky coordinates and the phase differences between the two orthogonal receptors give rise to polrotation.

CHALLENGES AND LIMITATIONS
A general mathematical foundation of polarization calibration has been provided by Hamaker, J. P. (2000) and a more recent review is available in Smirnov (2011). While a general prescription has been available, the complexity of the problem and its computation-heavy nature have restricted most commonly available implementations to make some simplifying assumptions. This is an active area of research, especially in view of the upcoming ambitious facilities like the SKA, and new algorithms and implementations are being developed by multiple groups across the world (e.g. Smirnov 2011;Mitchell et al. 2008;Salvini & Wijnholds 2014;Kenyon et al. 2018, etc.). This section lists the challenges of lowfrequency solar observation and the limitation of conventional algorithms and earlier attempts.

Challenges of Low-frequency Solar Observations
These challenges are related to the large FoVs and the nature of low-frequency aperture array instruments. Conventional methods correct for polconversion using observations of a strong unpolarized calibrator source. Correcting for polrotation requires observations of a strong source with known polarization properties (Hamaker, J. P. 2000;Hales 2017).
The large FoV of low radio frequency aperture array instruments imply that: 1. As there are no moving parts and the beam is steered electronically, the primary beam of instrument can vary dramatically with pointing direction.
2. Given the large FoV and nature of low radio frequency sky, there is no single polarized source strong enough to dominate the observed V ij .
This imposes the requirement that, for good polconversion calibration, the target field and the calibrator should be observed with the same pointing, which is rarely the case. The lack of a dominant polarized source makes it hard to do polrotation calibration. In addition, the low radio frequencies of observation require one to contend with the direction-dependent ionospheric distortion and Faraday rotation (FR). A comprehensive discussion of a successful approach to deal with these challenges has been provided by Lenc et al. (2017). In addition to the ones stated above, solar observations at these frequencies have additional challenges to deal with. These include: 1. The large flux density of the Sun and the wide FoVs imply that daytime calibrator observations are corrupted by solar contributions. Hence calibrators for solar observations are usually observed before sunrise or after sunset. The large time gap between the calibrator and the source observations, in addition to the difference in the pointing direction, reduce our ability to constrain the true state of the instrument and the ionosphere. This increases our reliance on self-calibrationbased methods to obtain high-dynamic-range solar radio images.
2. Solar emission is not assured to be unpolarized and can vary significantly across time and frequency. In the conventional approach, this makes the Sun unsuitable as a calibrator source for polconversion and polrotation.

Limitations of Conventional Algorithms
Most of the standard interferometric calibration and imaging packages like CASA and AIPS implement a linearized measurement equation. While they have been spectacularly successful in delivering high-quality polarimetric image, the following two assumptions must be satisfied for this formalism to be valid: 1. The instrumental polarization must be small (< 10%).
2. The fractional polarization of the sources used for calibration should be low (< 10%).
In the usual case of steerable antennas where the FoVs are small enough that one is never too far from the optical axis of the dish, the instrumental polarization usually meets this threshold. Also, the fractional polarization of standard polarization calibrators (e.g. 3C 286, 3C 138, 3C 48) is ≤ 10%. A linearized formalism is, hence, quite adequate for most applications. For our particular application of looking at the Sun with an aperture array, however, neither of these assumptions hold true. The instrumental polarization is a strong function of the pointing direction. It increases as one goes farther from the zenith and/or cardinal directions and is often much larger than 10%. The solar emission can vary dramatically in its intrinsic polarization from being unpolarized to being nearly 100% polarized. As calibrator observations are not possible during the solar observation with the MWA, one has to rely on the self-calibration using the Sun itself. The potentially very large fractional polarization of solar emission implies that one cannot rely on using them for polarization self-calibration with the linearized algorithms. This forces us to develop a more general formalism for polarimetric calibration. The linear approximation also restricts the dynamic range of the Stokes images (Smirnov 2011). In addition, in the linearized formulation, the instrumental leakages are calibrated only using the cross-hand visibilities (XY † , YX † or RL † , LR † ) (Hales 2017), and the parallel-hand visibilities (XX † , YY † or RR † , LL † ) are simply ignored. Hence, while the cross-hand visibilities are updated during polarization self-calibration, the parallel-hand correlations remain unchanged. As a consequence this algorithm is unsuitable for iterative implementation.

Previous Attempts and Their Limitations
While Mondal et al. (2019) have used self-calibrationbased methods with remarkable success to obtain highdynamic-range solar Stokes I images, their algorithm does not include polarimetric imaging. It also does not include absolute flux density calibration, which needs to be done independently (Kansabanik et al. 2022). Mc-Cauley et al. (2019) have demonstrated polarimetric so-lar radio imaging with the MWA. They used nighttime calibrator observations to estimate the instrumental and ionospheric gains. They used an ad-hoc approach to mitigate instrumental polarization, which we refer to as Method I in the following text. The assumptions and requirements of Method I, and the limitations they impose are listed below: 1. The leakage from Stokes I to other Stokes components remains essentially constant across the angular span of the solar disk. While this assumption is valid for some pointing, it is not true in general and this is demonstrated in detail in Sec. 4.5.2.
2. S Q = S U = 0, where S Q and S U represent the Stokes Q and U components of the solar flux density, i.e. the linearly polarized emission from the Sun is assumed to be exactly zero. While no linearly polarized emission has been reported from the Sun at low radio frequencies yet, the new generation of instruments can now provide spectroscopic snapshot images over spectral spans as small a few kHz, as opposed to order MHz available earlier. Assuming the linearly polarized flux density to be zero precludes their discovery and locks us out of an interesting discovery phase space.
3. Method I relies on the fact that the fractional circular polarization from the quiet Sun is expected to be small. It attempts to estimate an epochdependent instrumental leakage by minimizing the total number of pixels that show a fractional circular polarization larger than a chosen threshold, r c . It inherently assumes that this emission is coming from quiet-Sun regions. At low radio frequencies, the presence of multiple simultaneous active sources on the Sun can, however, limit the regions of the solar disk with quiet Sun emission (Mohan & Oberoi 2017). Even though the area occupied by the active regions is a small fraction of the solar disk at optical and higher frequencies, at low radio frequencies even the smallest active region gets scatter broadened to a few arcmin (e.g., Kontar et al. 2017;Mohan 2021b). In addition, the intrinsic brightness temperatures, T B , of various emissions often found to appear simultaneously on the Sun vary from 10 4 K for gyrosynchrotron emission to up to 10 13 K for type III solar radio bursts. The presence of very bright nonthermal emission has two consequences -one they lead to an increase in the system temperature and consequently the thermal noise in the image; and two they impose a larger imaging dynamic range requirement to be able to image the 10 6 K quiet-Sun regions in the presence of much brighter nonthermal emission.
The primary merit of Method I is that it enabled the authors to get to interesting science using an approximate and quick correction of instrumental polarization and circumventing the effort and complexities of developing and implementing a formally correct polarimetric calibration algorithm (McCauley et al. 2019;Rahman et al. 2020). The ionospheric phases during the solar observations are expected to be significantly different compared to those determined from nighttime calibrator observation. Hence applying the gain solutions from the nighttime calibrator observations limits the images to much poorer fidelity than the intrinsic capability of the data. These images cannot be used for reliable measurements of low levels of circular polarization. Method I is also known to give rise to some spurious polarization for very bright solar radio bursts (Rahman et al. 2020).

ALGORITHM
This section describes a robust formal polarization calibration algorithm that overcomes the shortcomings mentioned in Sec. 3.3 and enables high-fidelity polarimetric imaging. It builds on three pillars; i) selfcalibration, ii) availability of a reliable instrumental beam model, and iii) some well-established properties of low-frequency solar radio emission. For low-radiofrequency solar observations with aperture arrays, it is not feasible to obtain calibrator observations at nearby times with the same primary beam pointing as used for target observations. This algorithm is, hence, designed to not require any calibrator observations. We refer to this algorithm as Polarimetry using Automated Imaging Routine for Compact Arrays for the Radio Sun (P-AIRCARS). A detailed description of the algorithm is presented here, and that of its implementation as a robust unsupervised pipeline will be the subject of a forthcoming paper (D. Kansabanik et al. 2022, in preparation).

Full Jones calibration algorithm
In conventional polarization calibration tools like CASA, all four observed visibilities between an antenna pair are written as separate equations in terms of instrumental gains and leakages. These equations are approximated up to first-order terms in leakages and solved separately to obtain the instrumental parameters (e.g. Hales 2017). In full Jones calibration, the Measurement Equation is solved as a 2 × 2 matrix equation. From the Eq. 2 we can write the coherence noise, S, as, where ||.|| F represents the Frobenius norm (Horn & Johnson 1985) of the matrix, V ij is the model visibility, V ij the observed visibility, and J i and J j represent Jones matrices for antennas i and j respectively. The instrumental Jones matrices are estimated by minimizing S. Minimization of S leads to the matrix generalization of the conventional scalar calibration, which has been used in several standard interferometric software packages like CASA (McMullin et al. 2007), AIPS (Wells 1985), flagcal (Prasad & Chengalur 2012), and classical antsol (Bhatnagar & Nityananda 2001), where Jones matrices were replaced by a single complex number. The full Jones calibration was introduced by Hamaker, J. P. (2000) and Mitchell et al. (2008) and was later optimized as StefCal (Salvini & Wijnholds 2014). We have used the recently developed full Jones calibration software package, CubiCal (Kenyon et al. 2018;Sob et al. 2019), which uses complex optimization and the Wirtinger derivative (Wirtinger 1927). In brief, minimization of S reduces to an analytical update rule of J i in terms of V ij , V ij and the Jones matrices of other antennas as, The Jones matrices of all the antenna elements are initialized as the identity matrix. J i s are estimated using Equation 5 in subsequent iterations until the absolute value of the changes in the Jones terms and S fall below a small positive number ( ∼ 10 −6 ) between two consecutive iterations. We find that solutions generally converge within ∼20 -30 iterations.

Self-calibration algorithm of P-AIRCARS
In radio interferometry, it is standard practice to write the instrumental Jones matrices as a chain of independent 2×2 matrices, each with its distinct physical origin and referred to as the Jones chain (Smirnov 2011). In P-AIRCARS, the net Jones matrix for the ith antenna is given by, where ν, t and l refer to the frequency and time of observation and direction (θ, φ) respectively. These individual terms in Eq. 6 for antenna i are: represents the frequency-independent time-variable instrumental gain. This term also includes the directiondependent gains arising due to propagation effects. At low radio frequencies, they correspond primarily to ionospheric phase.
represents the timeindependent instrumental bandpass response.
is the Jones matrix representing the phase difference between the two orthogonal receptors for the reference antenna, given by ψ(ν, t).
the error on the ideal instrumental primary beam model.
The Sun being by far the dominant source in the sky implies that we are in essentially a small FoV regime. This offers the advantage that the direction dependence due to the ionospheric effects included on G i (t, l) can be ignored. We have verified that the direction dependence of E i (ν, t, l) cannot be ignored over the angular span of the Sun, but that of the D i (ν, t, l) can be, as discussed in detail in Sec. 4.6. Henceforth in this work, we ignore their direction dependence and regard If we write the V ij in terms the of sky brightness matrix, B( l) using the van Cittert-Zernike theorem (Thompson et al. 2017), and neglecting the noise term, Eq. 3 can be written as , where l, m and n are the direction cosines of l; and u, v, and w are the components of the baseline vector in units of the wavelength. Using Eqs. 6 and 7, we can write V ij as, Here the quantities outside the square bracket are the direction-independent Jones terms and the quantities inside the square brackets are the direction-dependent terms. We estimate each of these Jones terms step-bystep. The flowchart of the P-AIRCARS algorithm for estimating these Jones terms is shown in Fig. 1.

Intensity self-calibration
We first estimate the time-variable instrumental and ionospheric gain, G i (t), normalized over all antenna elements of the array. We choose a single frequency (ν = ν 0 ) channel for intensity self-calibration and set B i (ν 0 ) = 1. We write Eq. 8 as is the apparent model visibility for a single spectral channel after intensity self-calibration. Our intensity self-calibration algorithm follows the same philosophy as AIRCARS (Mondal et al. 2019), and our implementation incorporates additional improvements and optimizations. Making use of the compact and centrally condensed array configuration of the MWA and the very high flux density of the Sun, we estimate both antenna gains, G i (t) and V ij,Gcor (t) iteratively. At any given time stamp, t, we start with the calibration of only the phases of G i (t) using all baselines with at least one of their antennas from the dense central core of the MWA. Antennas at larger distances are progressively added. The process of intensity self-calibration with increasing number of distant antennas is similar to that implemented in AIRCARS (Mondal et al. 2019 Polarization self-calibration Figure 1. Flowchart describing the self-calibration algorithm of P-AIRCARS. The first major calibration block is intensity self-calibration, which is similar to that implemented in AIRCARS (marked in blue). The orange shaded box represents the polarization calibration block. The gray shaded box is a subset of the polarization calibration and shows the steps for polarization self-calibration. The Jones terms in Eq. 6 is being solved at each steps are denoted inside brackets.
success has already been demonstrated by many studies relying on high-dynamic range-imaging (e.g, Mohan et al. 2019b,a;Mondal et al. 2020b,a;Mohan 2021a,b;Kansabanik et al. 2022). The antennas forming the compact MWA core have the similar ionosphere along their line of sight (LoS) and thus have similar ionospheric phases. With the increasing distance from the core, the ionospheric phase for the antenna elements, with respect to a reference antenna in the core, grows. Since the MWA has a centrally condensed core with a large number of antenna elements and a rapid fall off in the density of antennas outside the core, the visibility distribution arising from the core is dominated by the baselines between the core antennas. This provides a certain level of coherency to the baselines originating from the core antennas. A detailed statistical demonstration of this and the advantages it brings is the subject of an independent publication (D. Kansabanik, 2022, in preparation).
The number of baselines originating from the core antennas is more than 70% of the all available baselines of the MWA. In a nutshell, the coherency of these large number of baselines makes the self-calibration process a well constrained problem, and allows the P-AIRCARS to start the intensity phase-only self-calibration process even without any calibrator source. At initial stages of the intensity self-calibration, the phases of the antenna gains of the farther away antennas are not that well-constrained, as only a small number of baselines involved contribute to the solutions. However, these initial gains are sufficiently close that they provide a good starting point for a self-calibration process. Once this process converges, baselines between increasingly distant antennas are included in small steps in subsequent self-calibration rounds. This is done to ensure that the problem is always well conditioned and most of the antenna solutions are already close to their optimal values. As more longer baselines participate in the calibration, it improves both the gain solutions for the antennas farther away from the core and the source model improve. This process is continued until all of the baselines are included and the final self-calibration iterations have converged. Once all antennas have been added and the phase-only self-calibration has converged, the phases of all antennas are deemed to be determined. Then the algorithm moves to amplitude and phase selfcalibration. Once amplitude-phase self-calibration using all the baselines has converged, gains for all the antennas are available. These gain solutions are applied and move to the next steps of the calibration process.
AIRCARS assumes V ij,Gcor (ν 0 , t)s to be unpolarized and effectively uses the same source model for the XX and YY polarizations. Unlike AIRCARS, no assumptions are made about the polarimetric properties of V ij,Gcor (ν 0 , t). A 2 × 2 matrix calibration is performed without any constraints on V ij,Gcor (ν 0 , t) except that G i (t) is assumed to be diagonal. Figure 2a shows the image after the first round of phase-only self-calibration, where the solar disk is rather distorted. The dynamic range of this image is only 24. The image after the first round of amplitude-phase selfcalibration is shown in Fig. 2b. The coherence of the array has improved remarkably, the solar disk is well formed, and the dynamic range has increased to 385. Figure 2c shows the final output of the intensity selfcalibration process. The improvement in image quality is self-evident, and the dynamic range has reached to 491. We note that prior to AIRCARS, the highest imaging dynamic range for solar imaging at meter wavelengths was a few hundred ( 300) and the imaging fidelity was too poor to be able to reliably detect features of strength few percent of the peak (Mercier & Chambe 2009). We have chosen a quiet featureless Sun for this illustration, the high-dynamic-range and highfidelity imaging of which continue to remain challenging at low radio frequencies even for the new generation instruments (Vocks, C. et al. 2018). The presence of solar activity makes the calibration and imaging problem easier, and dynamic ranges exceeding 10 5 have been achieved using AIRCARS (Mondal et al. 2019).

Bandpass self-calibration
As the instrumental bandpass amplitude and phase vary across frequency, bandpass calibration is required before combining multiple spectral channels to make an image. As AIRCARS was designed for spectroscopic imaging and the flux density calibration was done using an independent nonimaging technique , it did not need to include bandpass calibration. Conventionally, instrumental bandpass is determined using standard flux density calibrator sources with known spectra. Lack of availability of suitable calibrators pushes us to rely on bandpass self-calibration. This in turn required us to find a way to deal with the degeneracy between the instrumental bandpass shape and the intrinsic spectral structure of the source (Sun).
To avoid intrinsic spectral structure, we carefully choose sufficiently quiet times for the bandpass self-calibration from the initial flux-density-calibrated dynamic spectrum. We obtain the initial flux density calibrated dynamic spectrum using the nonimaging technique mentioned earlier, which is independent of instrumental gains and is computationally much faster. Though precise flux density calibration of MWA solar observations can be done using the method described in Kansabanik et al. (2022), it is only applicable to bandpass calibrated visibilities or images, which are not yet available at this stage in the calibration process. Bandpass self-calibration is performed for narrow bandwidths of 1.28 MHz at a time, referred to as a 'picket'. The spectrum of the quiet Sun can justifiably be assumed to be flat across a picket and apparent source visibility after bandpass calibration, V ij,Bcor (ν, t) = V ij,Bcor (ν 0 , t).
A single time slice (t = t 0 ) is chosen for bandpass self-calibration. Bandpass self-calibration starts with V ij,Gcor (ν 0 , t 0 ) as the initial model. Equation 10 over the band can then be rewritten as, is the apparent model visibility after bandpass self-calibration given as, We find that inter-picket bandpass phases for MWA can be modeled well by a straight line (Sokolowski et al. 2020). Inter-picket bandpass amplitudes show a more complicated variation and are corrected using an independent method described by Kansabanik et al. (2022). As expected, bandpass self-calibration improves the image quality, and the dynamic range increases from 431 ( Fig. 2c) to 793 (Fig. 2d).

Polarization self-calibration
This section describes the different parts of the polarization self-calibration algorithm marked by the orange shaded region in Fig. 1.

Cross-hand phase calibration
During the intensity and bandpass self-calibration, the phase of the reference antenna is set to zero for each of the orthogonal receptors. There is, however, an arbitrary phase difference between the two orthogonal receptors of the reference antenna. This cross-hand phase, ψ, is thus a single number and is applied to all antennas a b d c Figure 2. Improvements in image quality during the intensity and bandpass self-calibration. Images shown here are at 88 MHz with a frequency resolution of 160 kHz and temporal resolution of 2 s. The red contours levels are at 0.3, 0.6, 2, 8, 20, 40, 60, 80 % of the peak flux density. The black ellipse at the bottom left is the point spread function. a. Image after the first round of intensity self-calibration. Bright emission is present near the phase center, but the solar disk is distorted. b. Image of the first amplitude-phase intensity self-calibration. c. Image after the end of the amplitude-phase intensity self-calibration. The noise level decreases significantly in comparison. d. Image showing the effect of bandpass self-calibration. The image noise is reduced further after the bandpass self-calibration. The dynamic ranges of the images are 24, 385, 491, and 793, respectively. by a single Jones matrix K cross (ν, t). In case of linearly polarized receptors, K cross (ν, t) causes a leakage from Stokes U to Stokes V and vice versa.
There are, however, only a few bright polarized sources available for cross-hand phase calibration at low frequencies (Lenc et al. 2017;Lenc et al. 2018). In addition to the reasons given in Sec. 3.2, the requirement of a strong polarized point source at the phase center dominating the measured visibility, imply that the standard cross-hand phase calibration method available in CASA cannot be used for aperture arrays with large FoVs. Cross-hand phase cannot be determined using the selfcalibration-based approaches. We use an image-based method for calibration of ψ, tailored for low-frequency  aperture array instruments with large FoV (Bernardi et al. 2013). An observation of a linearly polarized source with sufficiently high rotation measure (RM) is chosen and the direction-independent instrumental gain and bandpass calibrations obtained from suitable unpolarized calibrator observation are applied. It is important to distinguish the true source polarization from the leakage from Stokes I due to the instrumental primary beam. The ideal primary beam corrections are hence applied to account for instrumental leakage. Even after the ideal primary beam correction, the residual leakages from Stokes I to other Stokes parameters still remain. For well-designed and well-modeled instruments, departures from nonorthogonality of the dipoles are small, which implies that the Stokes Q to Stokes V leakage is small. For the MWA it is small enough to be ignored (Bernardi et al. 2013;Lenc et al. 2017). When a linearly polarized emission passes through magnetized plasma, the polarization angle (χ) rotates. The observed χ ∝ λ 2 , where λ is the wavelength of observation. The proportionality constant is referred to as Rotation Measure (RM) and, like all propagation effects, is an integral along the entire LoS. It depends on the distributions of the LoS component of magnetic field strength and electron density. Multiple mediums (e.g., interstellar medium, interplanetary medium and ionosphere) with different RMs contribute to the total rotation. RM synthesis (Brentjens & de Bruyn 2005) is Fourier synthesis technique to separate out these RM components in the Fourier domain. The leakage flux from Stokes I to V can also be thought of as yet another medium contributing to the observed RM. This leakage flux must appear at the instrumental RM in Fourier domain, which is typically a few rad m −2 . As the linearly polarized source is chosen to be at an RM significantly higher than the instrumental RM, any Stokes V emission detected at source RM in the Fourier space can only arise due to leakages from components that rotate with the source RM. Hence the observed Stokes V emission at source RM must arise due to the leakages from Stokes U. Thus the Stokes U and Stokes V flux density are estimated from the image using RM synthesis. For linearly polarized receptors, K cross (ν, t) causes the leakage from Stokes U to Stokes V. We vary ψ between −180 • and +180 • and determine the value of ψ for which the spurious Stokes V emission is minimized.
It has been found that the ψ is extremely stable across both time and frequency for the MWA. Lenc et al. (2017) found that ψ essentially remains constant across the MWA frequency band from 80 − 300 MHz. The GaLactic and Extragalactic All-sky MWA (GLEAM) survey team (private communication with Xiang Zhang, ICRAR) have recently determined that ψ is stable over time scales of years. P-AIRCARS generally uses the values of ψ available a priori, but does provide the flexibility to estimate it from the nearest observations of a linearly polarized source and apply the corresponding correction. The extreme stability across time and frequency allows the use of nighttime observations from even months away to estimate ψ.
After the correction of the cross-hand phase, V ij,Xcor at any time, t and any frequency, ν, is obtained from Eq. 12 as

Ideal primary beam correction
For aperture array instruments, a significant part of the instrumental leakage comes from the pointingdependent instrumental primary beam. In reality, it is rarely feasible to measure the full-Stokes primary beam for aperture arrays. Hence, an ideal model of the primary beam, obtained using electromagnetic simulations is generally used. There are a few primary beam models available for the MWA -an analytical beam model (Ord et al. 2010); an average embedded element (AEE) beam model by Sutinjo et al. (2015); and full embedded element (FEE) beam models by Oberoi et al. (2017) and (Sokolowski et al. 2017). These models have steadily improved in their sophistication and performance. We have used the most recent ideal primary beam model for the MWA Sokolowski et al. (2017) for which the FEE beams have been simulated using the electromagnetic simulation tool, FEKO 2 . We note that P-AIRCARS allows the flexibility to use any primary beam model and will enable us to benefit from improved models as and when they become available. Naturally, the true primary beam response can deviate from the ideal beam model.
Method I assumes that the leakage from Stokes I to other Stokes components remains essentially constant across the angular span of the solar disk. While this assumption is valid for some pointings, it is not true in general. Percentage variation of this leakage across the solar disk is estimated to be as large as 50%. We define the percentage variation as, where l(θ, φ) is the leakage from Stokes I at the sky coordinate (θ, φ), l center is the leakage at the center of the solar disk, and ∆l(θ, φ) is the percentage change of l(θ, φ) with respect to l center . The variation of ∆l(θ, φ) over the Sun for two observing epochs is shown in Fig.  3. The Sun was close to half power point of the primary beam during the first epoch, 2014 May 04. During the second epoch, 2014 September 28, the Sun was close to the peak of the primary beam. We have found that l center is smaller for second epoch as compared to that for the first epoch. This is expected as the primary beam model is more accurate for near its peak. Nonetheless, the fractional variation of the leakages over the Sun for both the epochs are similar and not negligible. Hence, for precise polarization calibration, it essential to correct for the direction-dependent primary beam response. For an homogeneous array comprising identical antenna elements, the the modeled primary beam response can be assumed to be identical for all antenna elements. Thus we can substitute E i (ν, t,l) = E j (ν, t,l) = E(ν, t,l) in Eq. 13,  When a calibrator source is available, B(ν, t, l) is known and the only unknowns in Eq. 15 are the D i (ν, t) s. As no suitable calibrator observation is available, D i (ν, t) is not known a priori there is a degeneracy between B(ν, t, l) and D i (ν, t). A perturbative approach is used to break this degeneracy. As D i (ν, t) is small compared to E(ν, t,l)s, we approximate D i as identity matrix in Eq. 15 and obtain the source visibility, B 0 (ν, t, l), while incorporating the correction for E( l). Equation 15 is then takes the form, V ij,Xcor (ν, t) = E(ν, t,l) B 0 (ν, t,l) E † (ν, t,l) × e −2πi(uijl+vijm+wij(n−1)) dl dm n = B 0,app (ν, t, l) × e −2πi(uijl+vijm+wij(n−1)) dl dm n , where B 0,app (ν, t, l) is the apparent source visibility, and B 0 (ν, t, l) is the source visibility without the correction for D i (ν, t)s and E(ν, t,l). B 0 (ν, t, l) and B 0,app (ν, t, l) are related as B 0 (ν, t, l) differs from the true brightness matrix, B(ν, t, l), due to the following two reasons: 1. D i (ν, t) has been ignored in Eq. 15, which introduces errors in B 0 (ν, t, l).
2. In absence of calibrator observations or other independent astronomical constraints, the degeneracy between G i and V ij,Gcor cannot be broken unambiguously.
Although D i (ν, t) is expected to be small, this may not hold true for all pointings, especially the ones at low elevations. The top row of Fig. 4 shows the residual leakages after the ideal primary beam correction at 80 MHz for 2014 May 04. For this epoch the beam pointing was at the lowest permitted elevation for the MWA. The residual leakages for Stokes I to Stokes Q and U are ∼ 34 − 35% and Stokes I to Stokes V is ∼ 3 − 4%. For comparison, the bottom row of Fig. 4 shows the residual leakages for observation on 2014 September 28 at 119 MHz. The beam pointing for this epoch was close to the meridian, and the residual leakages are only a few percent. No systematic spatial variation of the residual leakages is seen in Fig. 4. The spatial variations over small angular scales are at 1-2% and are discussed in detail in Sec. 4.6. For precise polarization calibration, it is essential to reliably correct for D i (ν, t)s. The next section describes our approach for an image-based firstorder correction for D i (ν, t) s to B 0 (ν, t, l).

Image-based Leakage correction
As discussed in Sec. 2, the objective of polarimetric calibration is to correct for instrumental polconversion and polrotation, which cannot be achieved using self-calibration-based approaches alone. The most recent MWA primary beam model (Sokolowski et al. 2017) used in P-AIRCARS successfully reduces the instrumental leakages significantly. Though it can come close, no primary beam model can reproduce the true response exactly. To remove the residual instrumental leakages, we use some well-established physical properties of the quiet-Sun thermal emission at metee wavelengths to design an image-based correction. The characteristics of the quiet-Sun thermal emission we rely on are: 1. The brightness temperature of the quiet-Sun thermal emission is well known to lie in the range of 2. Quiet-Sun thermal emission has a very low level of circular polarization ( 1%), which arises due to propagation effects through the magnetized corona (Sastry 2009).

No linearly polarized emission is expected from the quiet-Sun thermal emission (Alissandrakis & Gary 2021).
No assumptions are made about the polarization properties of any active emissions as they depend on the emission mechanism, magnetic field strength, and topology, and are variable across time and frequency. Given these properties, one can argue that any Stokes Q and U emission seen in the quiet Sun region must arise due to residual instrumental leakages. Also, as the circular polarization of the quiet-Sun is 1%, the dominant contribution to the Stokes Q and U leakage must come from Stokes I. For the MWA, the misalignment of the dipoles with respect to the sky coordinates has been established to be small enough to give rise to insignificant mixing between Stokes Q and U (Lenc et al. 2017). An explicit correction for K cross (ν, t) has also been applied to correct for the mixing between Stokes U and V. Based on these arguments, instrumental leakages are corrected as follows: 1. A T B map is first made using the flux-densitycalibrated solar images (Kansabanik et al. 2022).
2. Regions where solar emission has reliably been detected are identified using an nσ lower threshold, where σ is the map RMS in a region far from the Sun, and n is usually chosen to lie between 10 and 6. This is the same threshold as used in the 'CLEAN'ing process during imaging and is determined and applied for each Stokes plane independently.
3. The regions lying between 10 5 and 10 6 K are considered to correspond to the quiet Sun.

Median values of Stokes Q and U fractions from
the quiet-Sun regions are computed and deemed to represent the leakages from Stokes I to Stokes Q and U, respectively.
5. The leakages thus determined are then subtracted from the Stokes Q and U maps of B 0 (ν, t, l).
We have only considered the leakages from Stokes I to Stokes Q and U. Stokes I to V leakages have generally been found to be consistent with 0% after the correction using the latest FEE beam model (Sokolowski et al. 2017). As discussed in Sec. 4.5.1 Stokes Q to Stokes V leakages are also small enough to be ignored for the MWA. Occasionally, however, when the beam pointing is at very low elevations, residual Stokes I to V leakage can grow to be as large as a few tens of percent. In such instances, an approach similar to what is used for estimating the first-order corrections for Stokes Q and U maps is used for obtaining a first-order estimate of Stokes I to V leakage. This is employed only when the median circular polarization in the quiet-Sun region is found to be > 2%. Method I used an older beam model (Sutinjo et al. 2015) which could not correct for the Stokes I to V leakage effectively, and needed to rely exclusively on this approach to estimate Stokes I to V leakage. In contrast, we use it sparingly, only when the Stokes I to V leakage is so large that our perturbative approach breaks down. These corrections account for polconversion, which arises due to the D i (ν, t) terms. There is little mixing between any of Stokes Q and U or Stokes Q and V for the MWA (Bernardi et al. 2013;Lenc et al. 2017). As K cross has also been corrected for, there is no mixing between Stokes U and V. Hence, no additional correction for polrotation is required. For aperture arrays, errors on the ideal beam are expected to be direction-dependent. This is true for the MWA as well (Lenc et al. 2017). We find that, although the leakages due to the ideal primary beam varies significantly over the angular extent of the Sun, the residual leakages due to D i (ν, t)s do not. Figure 4 shows the residual leakages from Stokes I to other Stokes parameters after the ideal primary beam correction. No significant systematic variations are seen across the solar disk after the subtraction of the median leakage, validating the assumption to treat the residual leakage or D i (ν, t)s as direction-independent quantities, as mentioned in Sec. 4.2. This image-based correction for leakages yields the source brightness matrix, B 1 (ν, t, l), which incorporates the first-order corrections for these leakages. Paralleling Eq. 16, B 1,app (ν, t, l) can be written as: We can now rewrite Eq. 15 using B 1,app ( l) and including the first-order correction of error matrices, D 1,i (ν, t)s, as V ij,Xcor (ν, t) = D 1,i (ν, t) B 1,app (ν, t, l) × e −2πi(uijl+vijm+wij(n−1)) dl dm n where V 1,ij,app (ν, t) is the apparent source visibility after first-order correction. D 1,i (ν, t)s can now be solved for iteratively using V 1,ij,app as the initial model, as discussed next.

Perturbative Correction : Residual Poldistortion Estimation and Correction
The V 1,ij,app (ν, t)s can be thought of as the full-Stokes sky model but with small errors arising largely from the deficiencies of first-order polarization calibration. The situation is analogous to the missing flux density problem in intensity self-calibration (Grobler et al. 2014), where it is addressed by using normalized solutions over the all antenna elements of the array. For the same reasons, a normalization factor also needs to be estimated in case of polarization self-calibration.
As the first-order corrections have already been applied in the process of determining V 1,ij,app (ν, t), D 1,i (ν, t) is expected to be small. D 1,i is estimated using full Jones matrix solver CubiCal. We define where D i (ν, t) represents the true values of the errors on the ideal instrumental primary beam mentioned in Eq. 15, and P D (ν, t) is the residual poldistortion left behind in the data after the first-order corrections. The mean of all D i (ν, t) is expected to lie close to the identity matrix and is the full-Stokes analog of the scalar normalization used in intensity self-calibration. For ease of notation we drop the explicit ν and t dependence of V 1,ij,app s, D 1,i s, D i s and P D s in the following text. We initiate the self-calibration process using V 1,ij,app s as the initial model. While D 1,i is not assured to be close to the Identity matrix, the first-order calibration already applied makes them small enough for a self-calibrationlike approach to converge. P D is estimated assuming that all D i are close to identity. To estimate P D , S D , the sum of the variance of D i with respect to identity matrix, I, is minimized. S D is defined as: For minimization, ∂SD ∂PD = 0 is imposed, leading to the following relation: We then correct each of the D 1,i for the P D and obtain D i as given by: Using Eq. 22, Eq. 19 can be written as, where V ij,app is the final apparent source visibility. The final true source brightness matrix, B(ν, t, l), is then obtained using the ideal primary beam as:

Estimation of Residual Leakages
For astronomical observations, residual leakages in the polarization images are usually determined using unpolarized celestial sources. The polarized emission detected from such sources, after polarization calibration, provides an estimate of the residual leakages from Stokes I to other Stokes parameters. For large FoV instruments, the leakage can be a strong function of direction. Despite this, an approach based on observations of multiple unpolarized sources has been successfully used (Lenc et al. 2017). They first determined the leakages toward the individual sources. Next a 2D polynomial of second order was fit to determine the leakages as a function of direction. It has recently been shown that numerous background sources can be detected in Stokes I with the MWA even in presence of the Sun (Kansabanik et al. 2022). However, to pursue a similar approach for determining leakages, a large number of these sources needs to be detected in other Stokes parameters as well, which is yet to be demonstrated.
We use the quiet-Sun regions in our images to estimate residual leakage fraction, L residual , as follows: 1. When Stokes Q, U and V emission is detected with more than 3σ significance over more than 50% of the quiet-Sun region, we use the median values of the pixels detected in each Stokes plane, and the define the residual leakages as: where med(L) Q , med(L) U and med(L) V are the median values of the Stokes Q, U, and V pixels detected with more than 3σ significance and I max is the maximum pixel value in Jy per beam in the quiet-Sun regions.
2. When no polarization is detected over the quiet-Sun region, we define a residual leakage limit based on the 3σ limit. In such situations we define the L residual as: where σ Q , σ U , σ V are the rms noise values in Jy per beam of the Stokes Q, U, and V images, respectively, close to the Sun and I max is defined as before.
3. Quiet-Sun emission in Stokes I may not always be detectable during the presence of very bright radio bursts (e.g., when a type III radio burst is in progress). In such cases, the residual leakage is estimated using the closest time stamp where the Stokes I quiet-Sun emission is detected.

A P-AIRCARS Stress Test
As a part of the process of development of P-AIRCARS and evaluating its efficacy, it was tested on some particularly challenging datasets. The most challenging of these observations comes from 2014 May 04 when the MWA was pointed to its lowest permissible elevation, where the sensitivity of the MWA and its polarization response are the poorest. This observation was at 80 MHz, where the flux density of the Sun is the lowest and Galactic background the strongest. Additionally, these data also correspond to quiet-Sun conditions,  where the Sun is essentially a featureless extended source and is hardest to calibrate and image. Processing these data, hence, correspond to a stress test of P-AIRCARS. It is quite reasonable to expect that, if P-AIRCARS can meet the challenge of calibrating and imaging these data, it will be able to successfully deal with most other MWA solar data. This section substantiates the performance of P-AIRCARS on these data. The Stokes images after the corrections of the modeled primary beam response are shown in the top row of Fig. 4. L residual for the Stokes Q, U, and V at this stage is ∼ 34, −35, and 2.8% respectively. As the Stokes I to V leakage was more than 2%, an imaged-based leakage correction was performed for Stokes V, as discussed in 4.6. The top row of Fig. 5 shows the percentage polarization over the quiet-Sun regions with more than 3σ detection in Stokes I after the image-based leakage corrections. The dynamic range of the Stokes I image is ∼500, and it is evident that Stokes I to Q, U leakages have reduced by about an order of magnitude (top panel of Fig. 5) beyond what was obtained by primary beam correction (top panel of Fig. 4). In this case Stokes Q emission is detected at > 3σ significance over the quiet-Sun region, but that at Stokes U and V is not detected. Hence we use Eq. 26 to obtain the limit of L residual for Stokes U and V. L residual ≈ −3.5 for Stokes Q calculated using Eq. 25 and L residual < |2| and |1|% for Stokes U and V, respectively, are calculated using Eq. 26. No systematic variation of the polarized emission is seen across the solar disk. After the image based leakage correction, several rounds of polarization self-calibration are performed and P D is corrected for at each iteration. This process is deemed converged when the absolute total polarized flux flux densities for Stokes Q, U, and V become stable. The rms of the Stokes images and the residual leakages are found to improve with every iteration of polarization self-calibration. The final Stokes Q, U, and V images are shown in the middle row of Fig.  5. Most of the regions in Stokes Q, U and V images are found to be consistent with noise. The L residual of the final Stokes Q, U, and V images have reached the values < |1.8|, |0.8| and |0.08|% respectively. The histogram of the pixel values of the Stokes Q, U and V leakage fraction are shown in the bottom row of Fig.  5. The median leakage values are already close to zero after image-based corrections, as expected, but tend to be asymmetric in some cases. In some cases, low-level artifacts are seen after image-based corrections. The regions with > 5σ detection after the image based correction are used for subsequent rounds of polarization self-calibration. The histograms of pixel values for the final images (shown in green) grow more symmetric and demonstrate a reduction in artifacts in the polarization images. The magnitude of this improvement can vary across different Stokes planes. For example, the histograms for Stokes Q, both after image-based correction (gray) and polarization self-calibration (green), are symmetric and very similar, demonstrating that imagebased correction was already very good and did not give rise to any significant artifacts. On the other hand, after image based correction, the Stokes V histogram shows a very skewed distribution with a positive tail extend-ing out to 7.5%. Polarization self-calibration leads to a much more symmetric distribution with a smaller span of ±2%. This demonstrates the efficacy of polarization self-calibration at reducing the artifacts in the Stokes V image. The situation for Stokes U is found to lie between these two regimes. While the improvements are always seen after polarization self-calibration, the trends seen here are not general and can differ from observation to observation.

P-AIRCARS uses the full Jones calibration package
CubiCal (Kenyon et al. 2018;Sob et al. 2019) for calibration and WSClean (Offringa et al. 2014) for imaging. The computational load of P-AIRCARS naturally depends on the details of the data, the choices made during analysis and computation hardware available. We present some numbers here to provide a general sense for the computational load and clock time taken for calibration. We use a 2 GHz CPU core with hyperthreading (two threads per core) as the benchmark device. The very first Stokes I self-calibration run on a chosen time and frequency slice (referred to as the reference time and frequency) takes about an hour. Once the gain solutions obtained from this reference slice have been applied to the dataset, the next step, bandpass self-calibration is performed in parallel on each of the 1.28 MHz wide 24 spectral chunks. This takes about 15-20 minutes per hyperthreaded core. The final step of polarization selfcalibration takes about 45 minutes per spectral slice and is also parallelized across the frequency axis. Thus for a typical P-AIRCARS run, full-Stokes calibration takes about 2 hours for a usual MWA dataset with 30.72 MHz bandwidth when parallelized across 25 cores (50 threads). As calibration (and imaging) are both done in a spectroscopic snapshot mode, only a tiny fraction of the entire dataset is needed for any individual calibration run making the memory footprint very small. More details about the implementation of the algorithm, the optimizations used, and the computational load will be presented in a forthcoming paper (D. Kansabanik et al., 2022, in preparation).

RESULTS AND DISCUSSION
The performance of P-AIRCARS on a particularly challenging quiet-Sun dataset has already been substantiated in Sec. 4.7.1. This section substantiates the various improvements that P-AIRCARS images represent and the science opportunities they enable. All of the solar data used here were taken under the MWA Project ID G0002. P-AIRCARS usually achieves < |2|% residual leakages for Stokes Q and Stokes U and < |1|% residual leakage for Stokes V. These values are comparable to a b Radio contour (80 MHz) Radio contour (80 MHz) Figure 6. Comparison between images made using AIRCARS and P-AIRCARS. The background images are LASCO white light coronagraph images. The red color map represents LASCO C2 images, and blue color map represents the LASCO C3 images. The green contours represents the radio images at 80 MHz. The contour levels are 0.2%, 0.4%, 0.6%, 0.8%, 2%, 4%, 8%, 20%, 40%, 60%, 80 % of the peak flux density. The filled ellipses at the lower left of the images are the PSF. a. The image made using AIRCARS has small extended emission from the CME as marked by the pink box. There are noise peaks at the 0.2% level. b. The image made using P-AIRCARS has extended emission over a larger region as marked by the yellow box. The radio emission covers the full white light CME.
what is generally achieved for high-quality astronomical observations (Lenc et al. 2017;Lenc et al. 2018;Riseley et al. 2018Riseley et al. , 2020 with the MWA and much better than those achieved by earlier spectropolarimetric solar studies (McCauley et al. 2019;Rahman et al. 2020). P-AIRCARS not only delivers a very small residual leakage, it also provides improved Stokes I imaging fidelity as compared to AIRCARS, as shown next.

Improvements in Stokes I imaging
In addition to the operational efficiency and polarimetric imaging capability, P-AIRCARS includes multiple improvements over the earlier Stokes I state-of-theart pipeline (AIRCARS; Mondal et al. 2019), which was focused on spectroscopic snapshot imaging. It is well unknown that the lack of calibration of instrumental polarization leakages leads to a reduction in the dynamic range (Bhatnagar & Nityananda 2001). Besides that, while imaging over multiple frequency channels, if the instrumental bandpass is not corrected, it leads to artifacts in the images. P-AIRCARS includes both these capabilities, which lead to a significant improvement in Stokes I image fidelity and noise properties.
A comparison of radio maps from AIRCARS (left panel) and P-AIRCARS (right panel) made using exactly the same data is shown in Fig. 6. The radio maps at 80 MHz have been superposed on Large An-gle and Spectrometric Coronagraph (LASCO Brueckner et al. 1995) images and use data of duration of 2 minutes and a bandwidth of 2 MHz. AIRCARS image shows a small weak emission feature from a CME on 2014 May 04 (marked by the pink box) over a region comparable to the point spread function (PSF).
The P-AIRCARS image shows emission at the similar strength but extended over a much larger region covering the full extent of the white light CME (marked by yellow box). It is evident that the imaging artifacts near the Sun in the P-AIRCARS image are at a lower level, and another weak extended emission feature lying on the western limb overlapping with the LASCO C3 FoV has a reliable detection. The noise characteristics of the P-AIRCARS image have also improved significantly. To substantiate this, Fig. 7 shows the entire FoV of the image shown in Fig. 6. It is evident that the noise characteristics of the AIRCARS image (left panel Fig. 7) is not uniform across the image, there is a bright and extended noise band running diagonally across the image from southwest to northeast. Dynamic range of this image is ∼ 1000. On the other hand, the noise characteristic of the P-AIRCARS image (right panel of Fig.  7) is quite uniform, and the large angular scale feature seen in AIRCARS image is no longer evident. The dynamic range of the P-AIRCARS image has improved by ∼ 80%, which leads to the value ∼ 1800. Despite its scientific merits being well recognized, gyrosynchrotron emission from CMEs has been only been detected in a handful of cases until date, of which an even smaller fraction are at meter wavelengths (Bastian et al. 2001;Bain et al. 2014). Using AIRCARS, Mondal et al. (2020b) have recently demonstrated the ability to detect gyrosynchrotron emission from CME plasma and model spatially resolved spectra to estimate plasma parameters of the CME for a slow and unremarkable CME. With further improved imaging from P-AIRCARS, it is now feasible to detect much fainter gyrosynchrotron emissions from CMEs and trace them out to higher coronal heights. This capability will make it possible to use this powerful tool for routinely measuring the CME plasma parameters and magnetic fields in the middle and upper corona.

Comparison with Method I
Section 3.3 discussed the earlier approach to spectropolarimetric solar imaging using the MWA, which has been referred to as the method I in this text (Mc-Cauley et al. 2019), and its limitations. The robust firstprinciples-based polarization calibration approach of P-AIRCARS, on the other hand, ensures that all known instrumental effects are corrected for. The Stokes images delivered by P-AIRCARS are limited primarily by the thermal noise of the data, and the rms noise seen in Stokes I and V images is significantly smaller than that seen in images from method I.
A comparison between images delivered by the method I (left panel) and P-AIRCARS (right panel) is shown in Fig. 8. Both images have been made from same observation on 2014 October 24, 03:46:00 UTC at 159 MHz. The red contours represent the Stokes I emission. The method I image shows artifacts and noise peaks at ∼2% of peak flux density. The largest noise peak in the P-AIRCARS image is at ∼0.8% of the peak flux density. The circular polarization percentage is shown by gray scales in regions where the Stokes I and V emissions are both detected with > 3σ significance. The dynamic ranges of both the Stokes I and V images have increased by an order of magnitude -for the Stokes I image from ∼ 35 to ∼ 577 and for the Stokes V image from ∼ 30 to ∼ 375.
This improvement in dynamic range enables the detection of a much weaker Stokes V emission with ∼3% circular polarization in the P-AIRCARS image (purple dotted box, right panel, Fig. 8), which was not de- The red contours represents the Stokes I emission. The contour levels are at 0.8%, 2%, 20%, 40%, 60%, and 80% of the peak flux density. Left panel: The image is made following the method I. The dynamic range of the Stokes I image is ∼35, and that of the Stokes V image is ∼30. There are noise peaks at the 2% level of the peak flux density in the Stokes I image. The residual Stokes V leakage is 5%. Right panel: The image is made using P-AIRCARS. The dynamic range of the Stokes I image is ∼577, and that of the Stokes V image is ∼375, which is an order of magnitude better than that for method I. The first noise peak appears at the 0.8% level of the peak flux density in the Stokes I image. We have detected another positive circularly polarized source (marked by the purple dashed box), which could not be detected in the image made using method I. This source has a circular polarization fraction ∼ 3%. The residual Stokes V leakage is 0.5%, which is also an order of magnitude improvement compared to that for method-I. tected using method I, and the region over which the circular polarization is detected with confidence in the P-AIRCARS image is substantially larger than that in the method I image. L residual has been estimated to be ∼ 5% for the method I and < |0.5|% for P-AIRCARS. The peak circular polarization fraction of the image from method I is ∼ −22% and from P-AIRCARS ∼ −27%. Considering the ∼ 5% uncertainty on circular polarization fraction of method I (McCauley et al. 2019), both these values are consistent. The very small L residual together with the high imaging dynamic range delivered by P-AIRCARS bodes very well for detection of the very low level of circularly polarized emission from the quiet-Sun thermal emission. This is also a very promising development for detection of polarized emission from gyrosynchrotron emission from the CME plasma.

Polarization Images of the Sun Using
P-AIRCARS Figure 9 shows the Stokes I and V images for example type I, type II, and type III solar radio bursts. These images are made at the native time and frequency resolution, 0.5 s and 40 kHz, of the observation. We have chosen different contour levels for different images to show the first noise peak of that image. Circular polarization fractions are shown in gray scale over the regions where both Stokes I and Stokes V emissions are detected with > 5σ significance. The peak circular polarization for these type I, type II and type III radio bursts are -89%, -53%, and -4.8%, respectively. The Stokes V flux densities of these bursts are ∼ 70, 180, and 70 SFU, respectively. In all cases we obtain images with DR varying between ∼ 500 and 1200. The high DR of these images enables us to detect a significant part of the quiet-Sun emission in Stokes I, even the presence of the bright active emissions. Interestingly, the peak of the Stokes V emission from the type III burst is slightly but clearly displaced toward the southeast from the peak of Stokes I emission, while they are coincident for the type I and type II bursts shown.
These, along with the quiet-Sun data showcased in 4.7.1 span a large range along the flux density axis. The disk-integrated flux density of the quiet Sun is ∼ 1.2 SFU, and peak surface brightness is 0.2 SFU per beam at 80 MHz. This demonstrates the capability of P-AIRCARS to produce high-DR images in a variety of different solar conditions. Additionally, the residual instrumental polarization leakages of 1% in all the Circularly polarized emission from a type I noise storm at 159 MHz. The dynamic range of the Stokes I and Stokes V images is 788 and 518, respectively. The red contour is at 0.8% of the peak Stokes I emission. The maximum circular polarization fraction is -89%. Middle panel: Circular polarization from a type II solar radio burst at 119 MHz. The dynamic range of the Stokes I and Stokes V images is 900 and 950, respectively. The red contour is at 0.5% of the peak Stokes I emission. The maximum circular polarization fraction is -53%. Right panel: Stokes V image of a type III solar radio burst at 118 MHz. The dynamic range of the Stokes I and Stokes V images are 1200 and 233, respectively. The red contour is at the 0.8% of the peak Stokes I emission. The maximum circular polarization fraction is -4.8%.
observations shown in Fig. 9, represent an improvement approaching an order of magnitude over the ∼ 5 − 10% leakages obtained earlier (McCauley et al. 2019;Rahman et al. 2020). To the best of our knowledge, these are the lowest to have been achieved for any meter-wavelength solar radio images and bring the exciting science target of measuring the weak circular polarization from the quiet Sun within reach.

Heliospheric Measurements Using Background Radio Sources
The Sun is the source with the highest flux density in the low-frequency sky, with even the quiet Sun spewing flux density exceeding 10 4 − 10 5 Jy in the MWA band. Apart from a handful of sources whose flux density goes up to hundreds of Jy, that of the bulk of other celestial sources lies in the range of a few Jy or weaker. Imaging a reasonably dense grid of background sources in the presence of the Sun in the FoV, hence, imposes a large dynamic range requirement, which had not been met until recently. The high dynamic-range of the images from P-AIRCARS now routinely allow us detect multiple background galactic and extragalactic radio sources even in the presence of the Sun. An example is shown in Fig. 7, which is made over 2.56 MHz and 10 s. The P-AIRCARS image shows a detection of 14 background sources with 5σ detection significance. The closest of these is at ∼20 R from the Sun and has a flux density of 4.9 Jy. Another example is available in Kansabanik et al. (2022), where the independently available flux densities of background sources were used for arriving at robust solar absolute flux density calibration for MWA solar observations. The rms noise in these images approaches that achieved by GLEAM, once the excess system temperature due to the Sun is taken into account (Kansabanik et al. 2022).
Measurements of interplanetary scintillation (IPS) of background radio sources have long been used to measure the electron density, properties of turbulence and velocities of CMEs, and solar wind in the heliosphere (e.g. Coles 1978;Manoharan & Ananthakrishnan 1990;Jackson et al. 1998) and more recently for driving the boundary conditions for magnetohydroynamic models of the solar wind (Yu et al. 2015). At the MWA, IPS observations have generally been done over small time windows while keeping the Sun at the null of the primary beam to avoid any contamination from solar emission (Morgan et al. 2018a,b;Chhetri et al. 2018). By providing the ability to routinely make Stokes I images of background sources, P-AIRCARS brings us a step closer to removing this limitation and opens the possibility of performing IPS observations without necessarily requiring to have the Sun in a null of the MWA primary beam.
While IPS is a very useful remote-sensing technique and provides information complementary to what is available from other observations, it is not sensitive to heliospheric magnetic fields: the key driver of space weather phenomena. By measuring the FR due to the heliospheric and/or CME plasma along the lines of sight to background linearly polarized sources or the diffuse galactic emission, radio observations provide the only known remote-sensing tool for measuring these magnetic fields. This approach has been successfully implemented by using radio beacons from satellites (e.g. Bird & Edenhofer 1990;Jensen et al. 2013;Wexler et al. 2019, etc.) as background sources, and more interestingly also using astronomical sources (e.g. Mancuso & Spangler 2000;Kooi et al. 2017Kooi et al. , 2021. These observation were carried out at higher frequencies and using small FoV instruments, which can sample only a small part of the heliosphere at any given time. Wide FoV instruments like the MWA can sample large swaths of the sky at any given time and can potentially track CMEs as they make way across the heliosphere. Measurements of FR simultaneously for a large numbers of pierce points across the CME/heliosphere, open the very exciting possibility of constraining the models for CME/heliospheric magnetic fields using these data Nakariakov et al. 2015). P-AIRCARS delivers precise polarization calibration and produces high-dynamic-range full-Stokes images and can already provide Stokes I maps of background sources. We are pursuing the objective of demonstrating making similar full-Stokes maps of background sources and enabling these heliospheric FR measurements.

CONCLUSIONS AND FUTURE WORK
We have developed a robust and comprehensive stateof-the-art polarization calibration algorithm tailored for the needs of low-frequency solar observations. P-AIRCARS builds on the learnings from our earlier Stokes I imaging pipeline (AIRCARS; Mondal et al. 2019) and uses the advantages endowed by the MWA design features to perform full polarimetric calibration without requiring dedicated observations of calibrator sources. The key MWA design advantages in this context are the dense and compact core of the MWA array layout and its simple and well-characterized hardware.
Together these ensure that nearby antennas, the ones looking through essentially the same ionospheric patch, maintain a good degree of coherence. A forthcoming paper will provide a detailed discussion and demonstration of these aspects (D. Kansabanik, 2022, in preparation). All the P-AIRCARS images shown here were made without using any calibrator observations. Polarization self-calibration was first demonstrated on simulated data by Hamaker (2006), but it had never been used for solar imaging. This work presents the first demonstration of solar polarization self-calibration and its ability to achieve high-dynamic-range and highfidelity full-Stokes images over a large range of solar conditions. The residual Stokes leakages for these images are on par with the usual astronomical images.
Though P-AIRCARS was developed with polarization calibration of the solar observations in mind, at its core the algorithm is general and does not impose any solar specific constraints. Its perturbative approach can be used for full Jones polarization self-calibration of the astronomical observations when a good initial sky model is available for a first-order calibration.
The perturbative algorithm used in P-AIRCARS works well for homogeneous arrays like the MWA, where the ideal primary beam response of all antenna elements are essentially identical. MWA antenna elements are made of a total 16 bow-tie dipoles arranged in a 4×4 grid (Tingay et al. 2013). The MWA beam is modeled assuming that all of the 16 dipoles are healthy (Sokolowski et al. 2017). It has been shown using satellite measurements that even when one or two of the dipoles fail, it does not change the primary beam response close to its peak in a significant manner (Line et al. 2018). However, for precise polarization calibration being pursued here, these small changes do need to be accounted for. Presently, in P-AIRCARS we reject all antenna elements with even a single dipole failure. Though it does lead to a loss of sensitivity, it is usually tolerable as the number of such elements is usually small. However for science applications close to the edge of the sensitivity limits, e.g. detection of CME gyrosynchrotron emission, it can become important to retain the sensitivity offered by the elements with defective dipoles. While the MWA beams can be modeled well for any subset of working dipoles (Sokolowski et al. 2017), it breaks the assumption of the array being a homogeneous one. The implication for P-AIRCARS is that an image-based approach for corrections for primary beams is no longer tenable (Eq. 16). One must then use a class of algorithms referred to in the literature as projection algorithms, which can correct for image plane effects in the visibility domain. These algorithms can be used for correcting artifacts arising from a wide range of causes, ranging from the so called w term to the antenna-to-antenna differences in primary beams even for an array with identical elements and ionospheric phase screens. The algorithm of relevance is the one referred to as the aw projection algorithm (Jagannathan et al. 2017(Jagannathan et al. , 2018Sekhar et al. 2021). It applies baseline-based corrections for primary beams in the visibility domain and is computationally very intensive. As efficient implementations of such algorithms become available and computational capacity available to us grows, it will become interesting to explore their use for scientifically interesting datasets with significant numbers of dipole failures to squeeze the most out of these data. P-AIRCARS has been developed with the future SKA in mind. It can easily be adapted for unsupervised generation of high-fidelity high-dynamic-range full-Stokes images for the SKA and other similar instruments with a dense central core. Our intent is to make a stable and mature implementation of P-AIRCARS available to the larger solar physics community. Producing high-quality solar radio interferometric images involves a steep learning curve, and its practise has remained limited to a small subset of the solar physics community. It is our belief that the lack of availability of a robust tool suitable the nonspecialist has long limited the use of radio observations in solar studies. We envisage that P-AIRCARS will prove to be a very useful tool for the solar and heliospheric physics community in times to come by filling this gap and making high quality full Stokes solar radio imaging accessible.