{'fraction_non_alphanumeric': 0.029595722160364, 'fraction_numerical': 0.010914539431179293, 'mean_word_length': 4.82636552939292, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Unified Calibration Framework for 21 cm Cosmology

Calibration precision is currently a limiting systematic in 21 cm cosmology experiments. While there are innumerable calibration approaches, most can be categorized as either `sky-based,' relying on an extremely accurate model of astronomical foreground emission, or `redundant,' requiring a precisely regular array with near-identical antenna response patterns. Both of these classes of calibration are inflexible to the realities of interferometric measurement. In practice, errors in the foreground model, antenna position offsets, and beam response inhomogeneities degrade calibration performance and contaminate the cosmological signal. Here we show that sky-based and redundant calibration can be unified into a highly general and physically motivated calibration framework based on a Bayesian statistical formalism. Our new framework includes sky and redundant calibration as special cases but can additionally support relaxing the rigid assumptions implicit in those approaches. Furthermore, we present novel calibration techniques such as redundant calibration for arrays with no redundant baselines, representing an alternative calibration method for imaging arrays such as the MWA Phase I. These new calibration approaches could mitigate systematics and reduce calibration error, thereby improving the precision of cosmological measurements.


INTRODUCTION
Measurement of the 21 cm cosmological power spectrum would constrain models of the Epoch of Reionization (EoR), Dark Ages, and Dark Energy. These measurements are contaminated by astrophysical foreground emission that is 4-5 orders-of-magnitude brighter than the cosmological signal. Separating the cosmological and foreground signals requires extremely precise instrumental calibration. As a result, development and characterization of precision calibration techniques for 21 cm cosmology has become an active area of research (Pen et al. 2009;Grobler et al. 2014;Newburgh et al. 2014; Barry et al. 2016;Berger et al. 2016;Grobler et al. 2016;Patil et al. 2016;van Weeren et al. 2016;Wijnholds et al. 2016;Ewall-Wice et al. 2016;Joseph et al. 2018;Li et al. 2018;Byrne et al. 2019;Orosz et al. 2019;Joseph et al. 2020).
In general, calibration approaches for cosmological 21 cm power spectrum measurements can be categorized as one of two types. 'Sky-based' calibration uses models of the sky and instrument to simulate data. Calibration then consists of fitting measurements to the simulation. In contrast, 'redundant calibration' relies on highly regular arrays with many redundant baselines measuring the same sky signal (Wieringa 1992;Parsons et al. 2010;Liu et al. 2010;Dillon & Parsons 2016;DeBoer et al. 2017;Dillon et al. 2018;Grobler et al. 2018;Kern et al. 2020;Dillon et al. 2020). Calibration then fits redundant measurements to each other, aiming for consistency between the measurements. Many diverse approaches fall into one of these two classes of calibration.
As the field of 21 cm cosmology pushes the limits of precision interferometric calibration it has become increasingly important that calibration frameworks mitigate error while capturing all instrumental systematics. Sky-based calibration approaches assume very good a priori models of the sky and instrument; model errors quickly degrade calibration and can preclude detection of the cosmological signal (Barry et al. 2016;Ewall-Wice et al. 2016). Redundant calibration is likewise vulnerable to sky model errors (Byrne et al. 2019) and experiences further calibration errors from instrument non-redundancies. Antenna position offsets and beam response inhomogeneities break redundant calibration's assumption of baseline redundancy and produce calibration errors (Joseph et al. 2018;Li et al. 2018;Orosz et al. 2019). The field requires novel calibration approaches that are resilient to sky and instrument model errors, can capture array non-redundancies, and mitigate contamination of the cosmological signal.
Combined calibration approaches integrate aspects of sky-based and redundant calibration. The Murchison Widefield Array (MWA) Phase II is a hybrid array that incorporates regular hexagonal sub-arrays within an otherwise pseudo-random array . Li et al. 2018 andLi et al. 2019 combine redundant calibration of the sub-arrays with sky-based calibration for the remainder of the array. Sievers 2017 proposes a calibration algorithm that relaxes redundancy requirements for redundant calibration of a regular array. This paper expands on these ideas to present a fully unified calibration framework that is statistically rigorous, highly flexible, and physically-motivated.
In §2 we introduce a calibration formalism based on a Bayesian statistical approach. In §3 we re-derive simple sky-based and redundant calibration frameworks with a focus on delineating the implicit assumptions of those approaches. In §4 we present a novel calibration framework that unifies redundant and sky-based calibration. Next, we explore extensions to this framework that allow for direction-dependent calibration ( §5), calibration across frequencies ( §6), and fully-polarized calibration ( §7).
Notation: Throughout this paper non-bold italicized variables represent scalars (e.g. f ), bold italicized variables represent vectors (e.g. v), and bold sans-serif variables represent matrices (e.g. A). Italicized subscripts represent indices (e.g. v j is the j-th element of the vector v) and non-italic subscripts are simply descriptive (e.g. u init ). For a summary of variables and expressions used throughout the paper see Table 1.

Bayesian Statement of the Calibration Problem
Calibration can be interpreted as a model fitting problem that fits tunable parameters to data. The fitting procedure asks, 'what are the most likely calibration parameters given the data?' If we define our visibility data as v and our tunable calibration parameters as θ, we can write the probability of θ given v as P (θ|v). This quantity is also called the 'likelihood function. ' We can use the likelihood function to calculate maximum likelihood parametersθ that maximize P (θ|v). As in most model fitting problems, there is no single 'right answer' for the form of the likelihood. One must construct a model that appropriately represents the data and mitigates systematics that contaminate the measurement. As precision calibration is crucial to the success of 21cm cosmology, choosing a well-motivated form for the likelihood function is extremely important. This paper describes a general calibration framework and explores different physically motivated likelihood functions that could improve calibration performance for 21cm cosmology.
From Bayes' theorem we can write the likelihood function as where P (v|θ) is the probability of the data v given model parameters θ. P (θ) is the prior probability, i.e. the probability of θ independent of the data. P (v) are is the marginal likelihood of v and is constant across all models. We can therefore simplify the problem by considering a proportionality P (θ|v) ∝ P (v|θ)P (θ) and maximizing the right hand side.
Maximizing the likelihood function P (θ|v) is equivalent to minimizing the negative log-likelihood. Calibration algorithms typical perform the latter procedure rather than explicitly maximizing P (θ|v). We describe the negative log-likelihood as where C 1 is an arbitrary positive constant and C 2 is an arbitrary constant of either sign. Neither constant affects the resultθ achieved by minimizing L(θ), and we therefore choose C 1 and C 2 to simplify the form of L(θ). For the case that P (θ|v) takes the form of a Gaussian distribution, the negative log-likelihood is equivalent to the least-squares cost function, also called the chi-squared and denoted χ 2 . For the sake of generality we do not assume that the likelihood is Gaussian. L(θ) can, but does not necessarily, take the form of a chi-squared. This formulation of the likelihood is completely general and θ can take any form. For example, traditional direction-independent calibration defines θ to be perantenna gains; direction-dependent calibration allows these gains to assume different values at different positions on the sky. Other calibration approaches such as redundant calibration may introduce additional tunable model parameters to parameterize uncertainties on the sky model. θ models the data according to some general function ζ(θ). Our calibration model assumes that this function can reconstruct the expectation value of the data for some values of the tunable parameters: where v is the expectation value of the data.
Variable or Expression Definition v data (for interferometric data, these are the measured visibilities) v expectation value of the data θ tunable calibration parameterŝ θ maximum-likelihood calibration parameters ζ(θ) function that maps the calibration parameters to the model of the data P (θ|v) likelihood function, equal to the probability that θ parameterizes the model given the data v P (v|θ) probability that the model parameterized by θ produces data v P (θ) prior probability of θ P (v) marginalized probability of v L(θ) negative log-likelihood; for a Gaussian likelihood function L(θ) is also called the chi-squared (or χ 2 ) I(θ) Fisher information matrix g antenna gains; these are an example of calibration parameters θ G(g) antenna gains written as a matrix; diagonal matrix with elements g a g * b u additional calibration parameters that are not g m the a priori estimator of u CT thermal covariance matrix of the data v σ 2 T elements of CT when CT is diagonal CM covariance matrix of the calibration parameters u σ 2 M elements of CM when CM is diagonal A matrix that maps calibration parameters u to the data x uv coordinate vector Bj(x) uv response of baseline j, assuming a continuous uv plane B matrix that encodes the uv responses of baselines f frequency γ a calibration parameters that parameterize the gain of antenna a across frequency; these are an example of calibration parameters θ V ab 2 × 2 visiblity matrix representing the 4 polarizations measured from baseline {a, b} J ak 2 × 2 Jones matrix for antenna a at sky location k S k 2 × 2 sky coherency matrix at location k E k polarized electric field vector of radiation emanating from sky location k η delay, Fourier dual of frequency f Once these choices are defined, calibration consists of fitting maximum-likelihood parametersθ. A final choice determines how these parameters are applied to the data to produce calibrated results. In some calibration frameworks only a subset of the fitted parameters get applied to the data. The rest of the fitted parameters are internal to the calibration optimization procedure. In general, calibration parameters related to the instrument response are involved in transforming uncalibrated data into calibrated data. Calibration parameters related to the sky model are internal to calibration.
All calibration frameworks described in this paper assume that the noise properties of the data are Gaussian. The probability distribution P (v|θ) from item (3) above therefore takes the form of a Gaussian distribution. While this is the most common choice, other calibration approaches relax the assumption of Gaussianity Ollier et al. 2017;Sob et al. 2020). The calibration frameworks in this paper can be easily extended to non-Gaussian distributions.

Quantifying Constraint of Calibration Solutions
A good parameterization of the model leads to highly constrained tunable parameters. The constraint on the tunable calibration parameters θ is quantified by the Fisher information.
The Fisher information matrix I is given by the negative Hessian of the log-likelihood: From Equation 3 we can also express the Fisher information in terms of the negative log-likelihood, L(θ): One can interpret the Fisher information as a measure of the curvature of the L(θ) hypersurface. At the maximum-likelihood point, where θ =θ, L(θ) experiences a minimum. We therefore expect that all elements of the Fisher information matrix are non-negative: Large values of this Fisher information indicate that L(θ) experiences a sharp minimum and the tunable parameters are highly constrained. Conversely, small Fisher information values mean the model is relatively agnostic to the parameter values and that the tunable parameters are not well-constrained. If I jk (θ) = 0 for all k then the parameter θ j is completely unconstrained by the model. It follows that the calibration solutions are degenerate. L(θ) does not experience a unique minimum but is rather minimized for any value of θ j . In this case, θ j is a degenerate parameter.
It is not always evident when calibration solutions are degenerate. Degenerate parameters can consist of linear combinations of the tunable parameters θ. We can calculate the number of degenerate parameters by taking the rank of the Fisher information matrix at the maximum-likelihood point. The calibration solutions are degenerate if I(θ) is singular. The number of degenerate parameters is equal to the dimensionality of the null space of I(θ); the degenerate parameters are the eigenvectors that span the null space.
Calibration degeneracies are a major challenge in 21cm cosmology. All degeneracies must be constrained to yield physical calibration solutions. Often, calibration consists of two stages where an initial calibration framework yields degenerate solutions and second step constrains those degeneracies (Liu et al. 2010;Zheng et al. 2014;Byrne et al. 2019;Kern et al. 2020).

TRADITIONAL DIRECTION-INDEPENDENT CALIBRATION APPROACHES
In this section we describe two simple directionindependent calibration frameworks, which we call 'skybased calibration' and 'redundant calibration.' We explain how these calibration frameworks connect with the statistical approach described in §2.
Both these calibration frameworks assume that the likelihood function is separable in frequency, time, and polarization. This means that each frequency channel, observation interval, and polarization mode can be calibrated independently, and we therefore omit explicit frequency, time, and polarization dependence. For a discussion of frequency calibration see §6; for a discussion of polarized calibration see §7. We assume per-time calibration throughout the paper.

Sky-Based Calibration
Traditional sky-based calibration makes the following choices: 1. It parameterizes the tunable instrument calibration parameters as a complex gain per antenna (and implicitly per polarization, frequency, and time step). θ = g where g has length equal to the number of antennas.
2. It models the data as where indices a and b index antennas. The combined indices ab indicate the index of the visibility formed by correlating signals from antennas a and b. m are model visibilities developed with a sky model and instrument simulator. The sky model typically consists of a point-source catalog, sometimes in conjunction with a diffuse foreground emission map.
3. It describes P (v|θ) as a Gaussian probability distribution: where the † symbol denotes the conjugate transpose. Here C T is the thermal covariance matrix, given by Furthermore, it assumes the visibilities are independent such that C T is diagonal. The probability function therefore takes the form where σ 2 T ,j are the diagonal elements of C T .
4. It uses a flat prior such that P (θ) = P (g) is a constant.
With these choices we can write the likelihood function for traditional sky-based calibration as Maximizing this likelihood function is equivalent to minimizing The resulting calibration solutions are degenerate. They have one degenerate parameter (implicitly per frequency, time, and polarization) that corresponds to the overall phase of the gains. One can see this degeneracy by noting that the transformation g → ge iφ does not change the form of L(g) for any real value φ. This phase can be interpreted as the absolute timing of incident radiation on the array and can be constrained by calibrating to a time-variable signal such as a pulsar.
All calibration frameworks described in this paper assume calibration to time-constant signals and therefore experience degeneracy in the overall phase of the gains. Typically this is constrained with a reference antenna. One then requires that Arg[ĝ ref ] = 0 whereĝ ref is the maximum-likelihood gain of the reference antenna and Arg denotes the complex phase. However, this means that systematic errors in the phase of the reference antenna affect all antennas. More sophisticated calibration approaches mitigate calibration errors in the reference antenna to better constrain the overall phase of the gains (Barry et al. 2019b;Li et al. 2019).
Traditional sky-based calibration assumes excellent knowledge of the sky and requires a highly complete sky model. It also assumes that the antenna responses are known up to a complex multiplicative factor. In other words, it assumes that the data can be fully modeled by the expression in Equation 8. This is clearly an inaccurate assumption. In reality, sky models are incomplete and inaccurate. These errors can propagate through the calibration process and degrade the calibration solutions (Grobler et al. 2014;Barry et al. 2016;Patil et al. 2016;Ewall-Wice et al. 2016;Joseph et al. 2020).
To account for sky model errors, we can replace the model visibilities with tunable calibration parameters u. Now θ = {g, u} and The likelihood function is now given by where P (u) is the prior on the tunable visibility parameters. Traditional sky-based calibration assumes that u is known exactly. Equation 15 converges to Equation 12 when where δ denotes the Dirac delta function. One expects that a more accurate calibration framework would not require so stringent a prior on u. In §4 we explore the implications of relaxing this assumption.

Redundant Calibration
We derive traditional redundant calibration (Wieringa 1992) from the following initial assumptions: 1. Like traditional sky-based calibration, traditional redundant calibration involves a tunable complex gain per antenna. However, redundant calibration assumes that the visibilities are not well-modeled and therefore also includes them as tunable sky model parameters. θ init = {g, u init } where u init corresponds to the visibilities. 2. We model the data as where a and b index antennas.
3. Like traditional sky-based calibration, traditional redundant calibration describes P (v|θ init ) as an independent Gaussian probability distribution as given by Equation 11.
4. As in traditional sky-based calibration, we use a flat prior for the gains g. The prior on the tunable calibration parameters therefore takes the form P (θ init ) = P (g)P (u init ) ∝ P (u init ).
Unlike traditional sky-based calibration, redundant calibration imposes a prior that visibilities from redundant baselines are correlated. If we assume the visibilities are Gaussian-distributed, then their prior takes the form where m init = u init . C Minit is the visibility covariance matrix and is given by Traditional redundant calibration assumes that redundant baselines measure the same sky. This is encoded in the covariance matrix by requiring that cov[u init,j , u init,k ] = var[u init,j ] if baselines j and k belong to the same redundant set and cov[u init,j , u init,k ] = 0 otherwise.
This appears to be a straightforward formulation of traditional redundant calibration -but it doesn't actually work. One cannot write down the negative loglikelihood L(g, u init ) based on the above assumptions because C M,init is singular: C M −1 ,init simply doesn't exist. Instead, we must remap u init to a new basis u that corresponds to an invertable C M . We find that an appropriate basis is to define u as the set of independent visibilities. u now has length equal to the number of redundant baseline sets rather then the total number of baselines. (See §A for a fully general discussion of working with singular covariance matrices.) Under this redefinition of u we get that where a and b index antennas and j indexes redundant baseline sets. The combined index ab denotes the baseline formed by correlating data from antennas a and b.
Here A is a rectangular matrix that maps u to the full set of visibilities. It has a number of rows equal to the total number of visibilities and number of columns equal to the number of redundant baseline sets. Its elements are given by Because traditional redundant calibration assumes that visibilities are uncorrelated between redundant baseline sets, C M is diagonal under this definition of u. We can therefore write the prior of u as where σ 2 M ,j are the diagonal elements of C M . Using this prior, the likelihood function becomes and L(g, u) is given by As with traditional sky-based calibration, minimizing this L(g, u) yields calibration solutions that are degenerate in the overall phase of the gains.
Traditional redundant calibration makes the further assumption that there is no prior knowledge of the sky. It takes the limit of a flat prior on u and lets σ 2 M ,j → ∞. This takes the second term in L(g, u) to zero, so This assumption of no prior sky knowledge introduces additional calibration degeneracies. In addition to the overall phase degeneracy of the gains, the calibration solutions are degenerate in the overall gain amplitude and the gradient of the gains' complex phase across the array (Liu et al. 2010). These degeneracies mean that redundant calibration must be combined with 'absolute calibration' to yield physical calibration solutions, where absolute calibration sets the values of the degenerate parameters. Minimizing Equation 26 is sometimes called 'relative calibration' to distinguish it from true redundant calibration, which must include both relative and absolute calibration steps (Zheng et al. 2014;Byrne et al. 2019;Kern et al. 2020). Traditional redundant calibration requires nearperfect redundancy of baselines within a redundant set. This means that antenna positions must be very close to a perfect grid and that each antennas' response must be near-identical up to a multiplicative gain. This can be very difficult to achieve in practice, where position errors and beam response inhomogeneities degrade the redundancy of an array (Joseph et al. 2018;Li et al. 2018;Orosz et al. 2019). In §4.5- §4.8 we propose new calibration frameworks that relax the requirement that redundant baselines be perfectly redundant.
Furthermore, traditional redundant calibration assumes a sparse array, such that the covariance between redundant baseline sets is zero. This is an appropriate assumption for some redundant arrays such as the Donald C. Backer Precision Array for Probing the EoR (PAPER; Parsons et al. 2010;Ali et al. 2015) and the hexagonal sub-arrays in the MWA Phase II configuration ), but it not a good assumption for compact redundant arrays such as the Hydrogen EoR Array (HERA; DeBoer et al. 2017) and the Hydrogen Intensity Real-Time Analysis Experiment (HIRAX; Newburgh et al. 2016). Because these compact arrays have large antenna collecting areas and small separations between antennas, they can have high covariance between baselines of similar but different lengths. In §4.6- §4.8 we relax the assumption that baselines in different redundant baseline sets have zero covariance. We expect that these more realistic assumptions will yield better calibration results for physical arrays.
Redundant calibration can be an attractive alternative to sky-based calibration because it has reduced reliance on the sky model. There is a misconception that one can use redundant calibration to calibrate interferometric data without any prior knowledge of the skyafter all, Equation 26 doesn't require a sky model at all! However, this neglects the role of the sky model in absolute calibration. Although the promise of redundant calibration is to reduce the impact of sky model errors on calibration solutions, it is inaccurate to assume one can calibrate without any prior knowledge of the sky. Furthermore, it is possible that introducing a sky model into the relative calibration step of redundant calibration as a prior on the fit visibility values could improve calibration performance. In §4 we explore calibration approaches that do just that.

A GENERAL FRAMEWORK FOR DIRECTION-INDEPENDENT CALIBRATION
In this section we describe a novel and highly general calibration framework. This framework allows for extensions to traditional sky-based and redundant calibration that relax some of the non-physical assumptions of those approaches.
The calibration framework described here is inherently direction-independent. For a discussion of directiondependent calibration see §5 or reference Kazemi et al. 2011Patil et al. 2016;van Weeren et al. 2016;Tasse et al. 2018;Albert et al. 2020;and Mertens et al. 2020. Additionally, it assumes per-frequency, per-time, and per-polarization calibration. §6 extends this framework to frequency-dependent calibration and §7 discusses fully polarized calibration techniques.

General Framework
In its most general form, this calibration framework makes the following assumptions: 1. It involves tunable complex gains per antenna (and implicitly per polarization, frequency, and time) g. The gains are not direction-dependent, meaning that this calibration framework is inherently direction-independent. In addition, this calibration framework involves tunable sky parameters u. The specific form of u depends on the class of calibration used. u can be interpreted as visibilities, pixels in the uv plane, source flux densities, or something else altogether. In general, u represents aspects of the sky model, possibly together with aspects of the instrument model, that are fit in calibration. The tunable calibration parameters are θ = {g, u}.
2. It models the data as Here A is a matrix that maps the sky parameters u to model visibilities. In more compact matrix multiplication notation, we can write where we have defined G(g) as a diagonal matrix with elements g a g * b . In the case that u represents visibilities, A is simply the identity matrix: A = 1.

It describes P (v|θ) as an independent Gaussian
probability distribution of the form of Equation 11.
4. It uses a Gaussian prior on u such that where m = u is a priori estimate of u. C M = cov[u, u † ] is an invertable matrix that encodes the covariances between elements of u. Furthermore, this calibration framework uses a flat prior on g and assumes that g and u are independent: Under this general framework, the likelihood function is given by and L(g, u) is given by C T is diagonal because the thermal noise on the visibilities is assumed to be independent. The calibration framework described here is highly abstracted. It represents a generalized treatment can be applied to practical calibration problems. In §4.2- §4.8 we delineate examples of calibration approaches that emerge from this framework. These are just some of the possible avenues for future exploration.

Sky-Based Calibration with Partial Sky Model Knowledge
Traditional sky-based calibration, presented in §3.1, assumes that the model of the visibilities is near-perfect up to a multiplicative complex gain. Of course this is not true. Model visibilities are susceptible to errors from incomplete knowledge of the sky and instrument.
We can make sky-based calibration more resilient to sky model errors by including tunable visibility parameters u in calibration. We assume that u is Gaussian distributed, so its prior is given by Equation 29. If we further assume that different visibilities are independent then C M is diagonal and the prior on u is where σ 2 M ,j are the diagonal elements of C M . We therefore get L(g, u) of the form and fit for calibration parameters θ = {g, u}. Note that this converges to traditional sky calibration (Equation 13) in the limit that σ 2 M ,j → 0. One common sky-based calibration technique is to calibrate only on selected baselines where the sky model is trusted. For example, this could involve calibrating only on baselines longer than 50 wavelengths in order to eliminate the short baselines that are poorly modeled by a point source catalog. We can replicate this technique from the formalism described in Equation 34 by taking the limits We now find that u j → m j for those baselines that are included in calibration: the fitted visibilities for those baselines are constrained to match the model visibilities.
On the other hand, u j for baselines that are excluded from calibration are completely unconstrained. They will therefore take on values such that g a g * b u ab = v ab for any values of the gains g a and g b .
While this calibration framework can replicate binary baseline selection it need not completely include or exclude baselines. σ 2 M quantifies the uncertainty on the model visibilities and can take any values. For example, it can be adjusted to be a function of baseline length to represent different levels of model confidence on different angular scales. Rather than completely eliminating short baselines from calibration, one could instead selectively downweight them, gradually increasing the value of σ 2 M ,j on subsequently shorter and shorter baselines. In fact, one can calculate σ 2 M empirically by measuring the agreement of the data and sky model as a function of baseline length. For example, the solid line in Figure  1 shows the power spectrum of data from a single twominute observation taken with the MWA Phase I. The data were processed through the Fast Holographic Deconvolution (FHD) 1 and Error Propagated Power Spectrum with Interleaved Observed Noise ( ppsilon) 2 software pipeline Jacobs et al. 2016;Barry et al. 2019a). The dashed line show the associated power spectrum of visibilities simulated with FHD from the GaLactic and Extragalactic All-sky Murchison Widefield Array (GLEAM) point source catalog (Hurley-Walker et al. 2017). The discrepancies between the data and simulation show the fidelity of the sky model. We see that the sky model is highly discrepant on large spatial scales while it is in good agreement with the data on small scales. We can set σ 2 M ,j proportional to the fractional difference between the data and sky model power spectra at the mode measured by baselines of length equal to that of baseline j. In this way we incorporate our knowledge of the sky model into our calibration model. We preferentially fit calibration solutions to the well-modeled long baseline measurements, but we also incorporate information from short baseline measurements to the extent that is appropriate for the given sky model.
In this section we have assumed independent visibilities: we do not account for covariances between visibilities. As we see with redundant calibration ( §3.2), this assumption does not hold for regular or close-packed arrays. In §4.4- §4.8 we explore calibration extensions that incorporate visibility covariances. Some calibration approaches in the literature consider the case of a point source sky model in which source intensities are not well-constrained (Mitchell et al. 2008;Sievers 2017). In this case, it could be advantageous to redefine u to be the source intensities rather than visibilities. Assuming that intensities are independent for each source, we can define L(g, u) to take the form

Sky-Based Calibration with Uncertain Source Intensitites
where m j = u j is the expected flux density of source j and σ 2 M ,j quantifies the uncertainty of that intensity. The matrix A models the instrument response. The element A jk is equal to the contribution to visibility j from a source of unity intensity at the position of source k.
Precision calibration for 21cm cosmology often includes tens of thousands of calibrator sources in the model. Allowing the intensities of all sources to vary is computationally infeasible. Instead, a more realistic approach would be to take σ 2 M ,j → 0, such that u j → m j , for the majority of sources j. Then σ 2 M ,j > 0 only for a select subset of troublesome sources that one expects are inaccurately modeled. This calibration approach is similar to directiondependent calibration ( §5). It can account for mismodeled apparent source intensities that vary across the sky. However, we classify this approach as directionindependent because the fitted gains -the calibration parameters that are actually applied to the data -are not direction-dependent. Importantly, this calibration framework does not allow different antennas to experience different intensity variations across the sky.

Redundant Calibration With a Sky Model and
Calibration of Hybrid (Redundant and Non-Redundant) Arrays As explained in §3.2, traditional redundant calibration does not require a sky model or an explicit instrument model in the relative calibration step. However, this means that it yields degenerate solutions. These degeneracies must be broken through an absolute calibration step that fits the degenerate parameters to model visibilities.
An extension to traditional redundant calibration incorporates model visibilities in L(g, u). This eliminates the need for a separate absolute calibration step and instead calculates optimal absolute calibration parameters in situ. It can also constrain the tunable visibility parameters u to approximate a known sky model.
From Equation 32 we get L(g, u) of the form As in §3.2 we define u to be the independent visibilities, with length equal to the number of redundant baseline sets. A maps the independent visibilites to the full set of visibilities and is given by Equation 22. Under those definitions C M is diagonal, so we can rewrite L(g, u) in the form of Equation 25: Here j indexes over redundant baseline sets. This calibration framework unifies the 'absolute' and 'relative' steps of redundant calibration. Like traditional sky-based calibration ( §3.1), it is degenerate only in the overall phase of the gains. Yet, like traditional redundant calibration, it incorporates information from baseline redundancy.
This framework also has interesting implications for 'hybrid' arrays, i.e. arrays with some redundant and some non-redundant elements. Phase II of the MWA, for example, contains two hexagonal sub-arrays that support redundant calibration. These redundant subarrays are contained within a larger pseudo-random array. Thus far, calibrating a hybrid array has required separate calibration operations for the redundant and non-redundant elements (Li et al. 2018(Li et al. , 2019. In order to exploit the redundancy of the sub-arrays one had to relatively calibrate those antennas with traditional redundant calibration techniques, therefore using baseline measurements only from within the sub-arrays. The calibration framework described in this section allows for calibration of the redundant sub-array antennas using all baseline measurements. This incorporates information from all baselines involving that antenna and can increase the calibration signal-to-noise.
To explore the scope of this calibration approach we describe two opposing limits of the calibration model. The first limit takes σ 2 M σ 2 T . In this limit we assume that the model visibilities m are well-known and require that u ≈ m. Extending this limit such that σ 2 M → 0 recovers traditional sky-based calibration.
In the opposing limit σ 2 M σ 2 T , where σ 2 T are elements of the diagonal matrix C T . The model visibilities m will then have a negligible contribution to the relative calibration parameters, i.e. the calibration parameters constrained by minimizing Equation 26. Instead, the model visibilities will fit the absolute calibration parameters only. This limit corresponds to a unified redundant calibration framework where absolute and relative calibration are encompassed in a single minimization problem. We note that extending this limit such that σ 2 M → ∞ recovers the relative calibration step of traditional redundant calibration.
In practice, neither of these limits are physicallymotivated. m cannot perfectly model the visibilities; on the other hand, it can give some information about the relative calibration solutions. The strength of this calibration approach is that it can inhabit the middle ground of some confidence in the model visibilities. As in §4.2, we can set σ 2 M empirically by measuring discrepancies between the data and model.

Redundant Calibration Accounting for Unmodeled Imperfect Redundancy
As explained in §3.2, we can describe traditional redundant calibration with a block-diagonal C M when the basis of sky parameters u corresponds to the full set of visibilities (denoted C M,init in §3.2). Traditional redundant calibration assumes cov[u j , u k ] = var[u j ] for visibilities u j and u k if baselines j and k belong to the same redundant set.
However, it is reasonable to assume that small antenna position errors and beam response inhomogeneities make cov[u j , u k ] < var[u j ] when j = k. We can account for this by applying a suppression term to the off-diagonal elements of C M . This suppression factor could be calculated empirically by measuring the covariance between visibilities from redundant baselines as the array measures different fields on the sky. It could also be baseline-or antenna-dependent. For example, if one antenna is known to have a particularly irregular beam response, the baseline covariance terms for baselines that include that antenna could be preferentially suppressed.
In general, suppressing the off-diagonal terms of a block-diagonal C M renders it invertable (in theory at least; in practice inverting C M could be computationally prohibitive). We may no longer need to recast u to a reduced basis such as the set of independent visibilities or the Singular Value Decomposition (SVD) basis described in §A. However, under this calibration approach the tunable visibility parameters u will be highly correlated. The maximum likelihood valuesû associated with redundant baselines will, by construction, have similar values. For this reason, we cannot use agreement ofû for redundant baselines as evidence of the baselines' true degree of redundancy.

Redundant Calibration with Modeled Imperfect
Redundancy §4.5 describes a calibration framework that accounts for unmodeled antenna position errors and beam response inhomogeneities. However, often antenna positions and beam responses can be measured to greater accuracy than they can be controlled. For example, an antenna in a redundant array may have a known position offset from its ideal position. In this case, the covariance matrix C M can be calculated from the modeled uv responses of the antennas. We represent the model of the uv response of baseline j as B j (x), where x is the uv position vector. These baseline response models can be developed from beam simulators or direct measurements and are readily available from instrument simulators such as FHD , pyuvsim 3 , OSKAR 4 , Precision Radio Interferometry Simulator (PRISim), or Common Astronomy Software Applications (CASA; see Jagannathan et al. 2017 for an example of using CASA with a fully-polarized primary beam model). The visibility u j relates to the true uv plane S(x) via From the definition of the covariance we get that If we assume no uncertainty on the uv response models, we can rewrite this as We now assume that different points in the uv plane are independent and have equal variance σ 2 M : We then get that This method of calculating elements of the covariance matrix allows one to accurately represent all amounts of baseline covariance in the calibration model. Redundant baselines have highly overlapping uv responses. In a closely packed array, baselines that are not constructed to be redundant may nonetheless overlap somewhat in uv coverage; this overlap introduces nonzero covariance in their measurements. Traditional redundant calibration ( §3.2) assumes that B j (x) = B k (x) when baselines j and k belong to the same redundant baseline set. It also assumes that when j and k do not belong to the same redundant baseline set their uv coverage does not overlap at all: If we instead calculate a covariance matrix from Equation 43 we can relax those assumptions and build a more physically-motivated calibration model that accounts for array non-redundancies and can incorporate covariances from baseline overlap in closely packed arrays such as HERA (DeBoer et al. 2017) and HIRAX (Newburgh et al. 2016).
In §4.2 we noted that a sky model may be more trusted for some baselines than others. The assumption in Equation 43 of uniform variance across the uv plane therefore must be relaxed. If we instead assume that the variance of the uv plane is approximately constant at scales equal to the size of a baseline response, we can rewrite Equation 43 as Here σ 2 M ,j is simply the variance of u j ; σ M,j is the standard deviation. As in §4.2, one could determine values of σ M empirically by comparing the data to the model, such as in Figure 1. Poorly-modeled power spectrum modes could be downweighted in calibration by increasing σ M,j for the associated baseline lengths. A baseline j could be excluded from calibration altogether by taking σ M,j → ∞.
Finally, one could combine the calibration approach described in this section with that of §4.5. One could expect that further non-redundancies exist beyond what is encoded in the baseline response model B j (x). In that case, one could represent non-modeled antenna position or beam response errors by suppressing the offdiagonal elements of C M . As in §4.5, this suppression factor could be estimated from empirical measurements of the covariance of visibilities. Calculating visibility covariances across observations of different parts of the sky could help validate or adjust a covariance matrix constructed from Equation 43.

'Redundant' Calibration of Compact Non-Redundant Arrays
The calibration formalism presented in §4.6 is, somewhat unexpectedly, fully generalizable to a compact non-redundant array. This means that one could implement 'redundant' calibration on measurements from an array that was not constructed to be redundant at all.
For example, Phase I of the MWA (Tingay et al. 2013) has no redundant baselines. As an imaging array, it was built to be pseudo-random with maximal uv coverage, in a sense making it as non-redundant as possible . At the same time, it has a highly compact core and is uv complete out to about 50 wavelengths. This high uv coverage means that one could construct a highly covariant C M from Equation 43 that would constrain visibilities across overlapping baselines.
Just as fully redundant arrays can be relatively calibrated without a sky model, as described in §3.2, compact arrays with very good uv coverage can be relatively calibrated without a sky model. One can calibrate these arrays by minimizing Equation 32 where C M is given by Equation 44. Taking the limit that σ 2 M σ 2 T reduces reliance on the sky model, instead constraining the relative calibration parameters from baseline covariance only. As in traditional redundant calibration, absolute calibration parameters would have to be constrained from the sky model. This could be a useful limit when the sky model is poorly constrained.
In this way the ideas behind redundant calibration can be extended to non-redundant compact arrays. We thereby reinterpret 'redundant' calibration as calibration that accounts for visibility covariances resulting from baseline's overlapping uv coverage. These calibration techniques have applications for many classes of interferometers, not just the highly redundant arrays typically associated with redundant calibration.

uv-Space Calibration
Traditional redundant calibration interprets the tunable sky parameters u to be visibilities, however one could instead reinterpret u as pixels of the uv plane. This is a natural parameter space for describing the model and means that C M does not depend on the instrument response. Instead, A encodes degridding by mapping the pixels of the uv plane to the visibilitites.
In §4.6 we describe the instrument response model as a continuous function in the uv plane, but in practice the uv plane is discretized. We can define an initial covariance matrix C M,init by rewriting Equation 43 for a discrete uv plane: where B is a rectangular matrix that maps uv pixels to visibilities. B has a number of columns equal to the number of pixels in the uv plane and a number of rows equal to the number of visibilities. The product BB † is equivalent to the holographic mapping function described in Sullivan et al. 2012. Now we define a new set of sky parameters u that correspond to pixels in the uv plane. B defines the mapping between u and u init : We now get that so the new covariance matrix is The fact that this covariance matrix is diagonal highlights the assumption of independent uv pixels in §4.6. Plugging this new parameterization of u into L(g, u) gives where m is the model of the uv plane and j indexes over uv plane pixels. We can also relax the assumption that σ 2 M ,j = σ 2 M for all uv pixels j and rewrite Equation 49 as Now σ 2 M can represent variations in the quality of the sky model across the uv plane. A point source sky model that omits diffuse emission could have σ 2 M increase for shorter baselines. Regions of uv space could be excluded from calibration altogether by taking the limit σ 2 M ,j → ∞ for those pixels j. As in §4.2, one can calculate σ 2 M empirically by comparing the model of the uv plane to measurements.
We call this new calibration framework 'uv -space calibration.' It offers a more elegant formulation of the calibration framework described in §4.6. Parameterizing u as uv pixels rather than visibilities highlights the statistical independence of the pixels by making C M diagonal. Here m depends on the sky model only and is calculated without an instrument simulator. The instrument model is moved from the second to the first term of L(g, u), where B encodes the baseline responses to the uv plane. Baseline covariances need not be calculated explicitly. Instead, this calibration framework implicitly constrains visiblities from fully-and partially-redundant baselines from their uv plane overlap.
uv -space calibration is a natural framework for unified calibration that models all baseline covariances. As with the calibration framework described in §4.6, it can account for small non-redundancies of a redundant array due to antenna position offsets and beam inhomogeneities. Like in §4.7, it models covariances stemming from baseline overlaps in the uv plane, thereby applying the ideas of redundant calibration to calibration of compact non-redundant arrays. Note that uv -space calibration incorporates modeled covariances only; it cannot include unmodeled imperfect redundancy as in §4.5. Finally, like the calibration framework described in §4.2, uv -space calibration can incorporate variable sky model uncertainties as a function of uv position to represent sky model incompleteness and preferentially fit calibration solutions to well-modeled uv modes.

DIRECTION-DEPENDENT CALIBRATION
Direction-independent calibration assumes that the shape of the antenna responses across the sky are wellmodeled a priori. In practice, unmodeled directiondependent effects can degrade interferometric performance. Direction-dependent calibration is a more flexi-ble alternative that can capture uncertainties in the spatial structure of the antenna responses.
Direction-dependent calibration has been explored extensively in the literature. In this section we show that generalized direction-dependent calibration can be formally described using the statistical framework presented in this paper. For a more thorough discussion of direction-dependent calibration see Kazemi et al. 2011Patil et al. 2016;van Weeren et al. 2016;Tasse et al. 2018;Albert et al. 2020;and Mertens et al. 2020.
A direction-dependent extension to traditional skybased calibration ( §3.1) uses the following assumptions: 1. Like all calibration frameworks presented in this paper, direction-dependent calibration uses tunable complex gains per antenna -but now, those gains can take different values at different positions on the sky. The tunable calibration parameters are given by θ = {g 1 , g 2 , . . . , g N } for N discrete sky sections, or facets. Here g n is the set of all gains for sky location n and has length equal to the number of antennas.
2. With direction-dependent calibration the visibility model ζ(θ) no longer takes the simple form of Equation 8. As in §4, we can represent the mapping from sky parameters to visibilities with a matrix A, but A is now a function of the directiondependent gains. We can write the visibility model as ζ(g 1 , g 2 , . . . , g N ) = A(g 1 , g 2 , . . . , g N )m. (51) Here m represents the sky model. It could be parameterized as point source intensities or the intensities of pixels on the sky. Unlike in traditional sky-based calibration, m cannot be parameterized as the estimate of the visibilities because the visibilities cannot be modeled a priori. Simulated visibilities will depend on the direction-dependent gains.
3. As in §3.1, one can describe P (v|θ) as an independent Gaussian probability distribution (Equation 11). This framework can be extended to other distributions with no loss of generality. 4. One could use a flat prior for the directiondependent gains. Alternatively, a common prior in direction-dependent calibration requires that gains from adjacent sky sections have similar values.
Direction-dependent calibration can be a good calibration approach when antenna responses are difficult to model. The direction-dependent gains can fit response shapes that are not well-modeled and vary antennato-antenna. Additionally, direction-dependent calibration can be important for modeling ionospheric effects. This is particularly critical for large arrays where antennas have different lines-of-sight through the ionosphere. Refraction through the ionosphere imposes directiondependent image distortions that vary in time as ionospheric conditions change (Jordan et al. 2017;Albert et al. 2020).
Direction-dependent, per-frequency calibration increases the calibration degrees-of-freedom in a way that can suppress diffuse foreground emission and the cosmological signal itself, leading to signal loss. One can protect against this signal loss by reducing calibration degrees-of-freedom across frequency (Patil et al. 2016). In §6 we explore alternatives to per-frequency calibration. With any calibration approach, one must quantify and mitigate potential signal loss stemming from coupling of the tunable calibration parameters and the cosmological signal.

FREQUENCY CALIBRATION
In the previous sections we implicitly assume a calibration likelihood function that is separable in frequency: (52) A likelihood function of this form allows for each frequency to be calibrated separately but precludes the calibration model from incorporating cross-frequency covariance information. In this section we describe calibration extensions that break the assumption of separability of the likelihood function in frequency. Frequency-dependent calibration errors pose a major problem for 21 cm cosmology. In particular, fast frequency structure in the calibrated gains can contaminate the cosmological signal by limiting the separability of the cosmological signal and the astronomical foregrounds.
Gain-based calibration assumes that the instrumental response is well-modeled as a multiplicative factor on the true sky signal. In other words, a visibility is given by v where u ab (f ) is the true sky signal at frequency f and G ab (f ) is the instrumental gain. If calibration is errorfree, one can divide the data by the gains to recover the true sky signal u ab (f ), which contains both the cosmological signal and contamination from foreground emission. The basis of 21 cm cosmology is the principle that these signals are spectrally distinguishable. Foreground contamination is assumed to be spectrally smooth and can therefore be separated from the specrally variant cosmological signal. If we use u ab (η) to represent the Fourier transform of the signal -where η is the Fourier dual of frequency, or delay, with units time -we expect u ab (η) to be free from foreground contamination at large η . The precise contaminated modes depend on the baseline length, and this contaminated region is often called the 'wedge' (Furlanetto et al. 2006;Datta et al. 2010;Morales & Wyithe 2010;Vedantham et al. 2012;Morales et al. 2012;Parsons et al. 2012;Trott et al. 2012;Dillon et al. 2013;Hazelton et al. 2013;Thyagarajan et al. 2013Thyagarajan et al. , 2015.
From the convolution theorem we can rewrite Equation 53 as where * denotes the convolution. Now consider the case that the gains calculated in calibration have an error that contributes to a particular delay mode η 0 . For that mode, G ab (η 0 ) → G ab (η 0 ) + ab (η 0 ), where ab (η 0 ) is the error on the gains. Now calibration cannot recover the true sky signal; instead, it recovers the signal ab (η) * u ab (η). The error on the gains therefore convolves all the true sky power, including the foreground power. If η 0 is large, this means that foreground power is coupled into the modes we expect to be free from contamination, otherwise known as the 'EoR window.' This contaminates the cosmological signal and can preclude a detection. For this reason, it is critical that calibration is essentially free from errors with fast frequency structure. These errors can take two forms. First of all, we require that calibration does not introduce new fast frequency structure by fitting thermal noise or sidelobe noise from errors in the calibration sky model. Secondly, we require that calibration accurately captures any true frequency structure in the instrumental response.

Parameterizing Gains Across Frequency
One approach to eliminating frequency-dependent calibration errors is to fit calibration solutions across frequency after-the-fact (Barry et al. 2019b;Li et al. 2019). With this approach, calibration consists of two steps. First, one calculates per-frequency calibration solutions. Next, the gains are adjusted to remove false frequency structure and fit the allowable degrees of freedom determined by a model of the instrument bandpass response.
A more optimal way of calibrating with a modeled bandpass is to fit the parameters of the bandpass directly (Yatawatta 2015;Mertens et al. 2020). We can model the gains as functions of tunable parameters γ a for each antenna a such that Here γ a could include the amplitudes and positions of known features in the bandpass or any other parameters related to fitting the instrument's bandpass response. If the bandpass is parameterized as a low-order polynomial, γ a could be the polynomial coefficients: where n max is the maximum mode in the low-order polynomial fit.
Because the parameters γ a are not per-frequency, calibration can no longer be parallelized across frequency. Instead, the algorithm must access data from across a large frequency range at once. This can prove computationally infeasible without novel calibration algorithms that allow for distributed optimization without sacrificing cross-frequency performance (Yatawatta 2015).
An extension to traditional sky-based calibration (Equation 13) that parameterizes the gains across frequencies has the following features: 1. It defines tunable calibration parameters corresponding to parameterizations of the antenna gains across frequency. θ = {γ 1 , γ 2 , . . . } where for antenna a the gain g a (f ) = g a (f, γ a ).
2. It models the data as where indices a and b index antennas.
3. As in previous sections, it describes P (v|θ) as an independent Gaussian probability distribution.
This version of sky-based calibration has a negative log-likelihood of the form Calibration consists of minimizing this quantity by varying the parameters γ a for each antenna a. This calibration approach allows for variable weighting across frequency. For example, frequency channels with greater noise contamination could be downweighted by increasing σ 2 T (f ) with respect to the uncontaminated channels. Frequency channels can be removed from calibration altogether by taking the limit σ 2 T (f ) → ∞. If the calibration parameters are chosen to span frequencies then they will interpolate the gains across the masked frequency channels.
This approach to fitting the gains is not only applicable to sky-based calibration. One can choose to parameterize the gains across frequency in any calibration framework.

Capturing Redundancy Across Frequency
Calibrating across frequencies allows redundant calibration to capture cross-frequency visibility covariances (Dillon & Parsons, private communication). It is possible that baselines of vastly different physical lengths are nonetheless highly covariant at different frequencies. A long baseline could measure the same sky modes at a low frequency that a short baseline measures at a higher frequency. This has interesting implications for redundant calibration.
To illustrate cross-frequency redundant calibration we consider an extension to redundant calibration as described in §3.2. Under this extension, calibration includes the following assumptions: 1. As in §3 and §4, it parameterizes the gains as per-antenna, per-frequency, and directionindependent.
(Alternatively, one could combine this calibration approach with the crossfrequency gain parameterizations described in §6.1. However, for simplicity we will describe a per-frequency gain parameterization.) It parameterizes the sky model with per-frequency visibilities u(f ). The tunable calibration parameters are θ = {g(f ), u(f )}.
2. As in §4, it models the data as for some matrix A that maps u to visibilities.
4. It uses a Gaussian prior on u and a flat prior on g such that C M now encodes covariances between every visibility at every frequency, with elements given by where j and k index baselines and f 1 and f 2 index frequency channels.
L(g, u) from Equation 32 now takes the form This converges to traditional redundant calibration in the limit that C M,jkf1f2 → 0 for f 1 = f 2 . Assuming a stable sky as a function of frequency, a physically-motivated construction of C M could follow the approach described in §4.6. From Equation 43 we get where B j (f, x) is the uv response of baseline j and S(f, x) is the uv plane at frequency f . If, as in §4.6, we assume that different points in the uv plane are independent and that the variance across the uv plane is constant, we get that Here σ 2 M,f f is the variance of uv pixels at frequency f ; σ 2 M,f1f2 is the covariance of uv pixels at frequencies f 1 and f 2 .
It follows from Equation 64 that the construction of a physically-motivated cross-frequency C M requires accurate modeling of the beam responses at all frequencies in question. In practice, beam responses vary considerably as a function of frequency. While baselines of different physical lengths may sample the same uv locations at different frequencies, we cannot expect their beam responses to be the same. Therefore, redundancy is degraded across frequencies. In other words, we expect the largest elements of C M,f1f2 to occur when f 1 = f 2 .
If we assume the sky is constant across our frequency range then σ 2 M,f1f2 = σ 2 M,f1f1 at all frequencies. However, in practice the sky has some frequency-dependence. We therefore expect that a more accurate model sets σ 2 M,f1f2 < σ 2 M,f1f1 when f 1 = f 2 . This further suppressed the elements of C M that correspond to covariances of visibilities at different frequencies.
Incorporating nonzero cross-frequency visibility covariances in calibration could be a powerful tool for constraining the frequency structure of the instrument response. However, beam shape variation across antennas and frequencies, along with frequency-dependence in the sky signal, suppress the covariance of visibilities at different frequencies. A realistic calibration model must account for these effects.

POLARIZED CALIBRATION
Polarized interferometry can be accomplished when each antenna has two polarization modes. This translates to four measured polarization modes from each baseline. The polarized measurement equation is given by where k indexes positions on the sky. Here V ab is a 2 × 2 matrix with elements corresponding to the four polarization modes measured by baseline {a, b}. The off-diagonal elements of V ab correspond to the 'crossvisibilities' formed by correlating measurements from different polarizations of the two antennas. J ak and J bk are Jones matrices associated with antennas a and b, respectively, and N ab gives the noise on the measurements. S k is the coherency matrix, equal to where E k is the electric field vector on the sky at position k Hamaker 2000). S k is diagonal for an unpolarized sky, but it does not follow that V ab will be diagonal, even for an instrument with orthogonal antenna polarizations. Although an antenna's two polarizations may measure orthogonal polarization modes of incident radiation from zenith (or, for a mechanically steerable antenna, from the pointing center), it will measure non-orthogonal polarization modes of off-zenith radiation. This effect is particularly pronounced for widefield instruments and fields with bright off-zenith sources. For this reason, one can expect that even unpolarized incident radiation couples into the measured cross-visibilities of an interferometer.
Equation 65 is often rewritten to describe the visibilities as a vector of length 4 rather than a 2 × 2 matrix. In that notation, Equation 65 becomes where J ak ⊗ J * bk is the 4 × 4 Kronecker product of the Jones matrices. v ab , s k , and n ab are the vector representations of the visibilities, the coherency matrix, and the noise, respectively.

Per-Polarization Calibration
Calibration frameworks presented thus far in this paper have assumed that L(θ) is separable in polarization. Under this assumption one can write the likelihood function as where p indexes the instrumental polarization modes.
Here θ p is the set of calibration parameters associated with instrumental polarization p. v ab,pp is the the (p, p) element of the matrix V ab ; v pp is the vector formed from those elements for all baselines. Note that crossvisibilties are excluded from per-polarization calibration, which amounts to leaving potentially valuable information on the table.
We can rewrite traditional sky-based calibration (Equation 13) with explicit polarization dependence as Here g a,p is the gain of the p-polarization of antenna a and m ab,pp is an estimate of the visibility v ab,pp . The model visibilities are produced from estimations of the sky coherency matrices S k and models of the Jones matrices J ak . Per-polarization calibration has a degeneracy associated with the average phase between polarizations. One can see this degeneracy by noting that the transformation g p → g p e iφ does not change the form of the L(g), where g p is the set of all p-polarized gains. This degeneracy means that per-polarization calibration cannot calibrate the cross-visibilities. One therefore cannot produce polarized images from per-polarization calibration.

Polarized Calibration with Per-Polarization Instrumental Gains
One approach to fully polarized calibration retains the per-polarization parameterization of the gains but incorporates the cross-visibilities. A polarized extension to traditional sky-based calibration would use L(g) of the form where p and q index instrumental polarizations. Because this calibration framework includes the crossvisibilities it does not have the degeneracy associated with an overall phase between polarizations. Instead, it has just one degenerate parameter corresponding to the overall phase of all gains. Although the gains are parameterized per-polarization, this calibration approach is fully polarized and can be used to produce polarized images.
Equation 70 represents a polarized extension to traditional sky-based calibration (Equation 13), which we have implicitly defined to be per-polarization. This polarized calibration approach could be combined with other calibration frameworks to achieve, for example, polarized redundant calibration (Kohn et al. 2019).

Calibrating Polarization Leakage
Per-polarization gain parameterization does not allow calibration to correct for polarization leakage. Polarization leakage occurs when power modeled in one visibility mode is measured as another . This can occur through a number of mechanisms. An antenna can be physically mis-aligned such that the so-called East-West polarized feed is somewhat skewed North-South. Antennas can experience cross-talk -for example, a signal from a North-South polarized feed could get coupled into the East-West feed's signal path. Additionally, the ionosphere can Faraday rotate one polarization mode into another. One can expand the calibration degrees-of-freedom to capture polarization leakage.
To correct for polarization leakage in calibration, we can parameterize each of the per-antenna gains as a 2 × 2 matrix with two polarization indices. g a,pr is then the gain of the p-polarization of antenna a with respect to incident r-polarized radiation. When p = r, g a,pr encodes the polarization leakage from mode r to mode p. An extension to traditional sky-based calibration that uses that calibration approach has L(g) of the form L(g) = ab pq rs 1 σ 2 T,abpqrs v ab,pq − g a,pr g * b,qs m ab,rs 2 .
(71) Here p and q index polarization modes of the measured visibilities and r and s index polarization modes of the modeled visibilities.

CONCLUSION
Precision calibration is essential to the success of 21 cm cosmology experiments. Low-level calibration errors contaminate the cosmological signal by limiting the separability of bright astrophysical foregrounds. As current experiments approach a detection of the 21 cm power spectrum from the EoR, novel calibration approaches must achieve new levels of precision.
While a diversity of interferometric calibration techniques exist, the current paradigm in the field delineates between sky-based and redundant calibration approaches. Sky-based calibration assumes very good a priori models of the sky and instrument response. It further assumes independent visibility measurements, neglecting visibility covariance from baseline response overlap in the uv plane. On the other hand, redundant calibration assumes that baselines within a redundant set are perfectly redundant both in antenna positions and response shapes, that baselines from different redundant baseline sets exhibit no covariance, and that a sky model can accurately constrain the absolute calibration parameters. The literature shows that violating these assumptions produces potentially catastrophic calibration errors. For example, Barry et al. 2016  The standard assumptions of sky-based and redundant calibration can be inflexible to demands of calibrating real data from physical interferometers. Nextgeneration calibration techniques, such as those presented in Sievers 2017 and this paper, combine elements of these disparate calibration approaches. New calibration frameworks must be statistically rigorous, physically-motivated, and adaptive to the complex physical systems they model. This paper presents a highly general statistical language that clarifies the assumptions implicit in a calibration formalism. It then describes new approaches to calibration that incorporate imperfect sky models, relax redundancy requirements, and apply the ideas of redundant calibration -namely, that baselines can be covariant with one another -to a wider class of arrays than those that are typically redundantly calibrated.
While the calibration frameworks described in this paper are highly general, further extensions to this work could expand upon them. One such extension is explicitly time-dependent calibration, where the system's evolution in time is parameterized in the calibration model. Another avenue of exploration could involve developing new approaches for constraining the overall phase of the calibration solutions. These could incorporate gated pulsar measurements, thereby augmenting skybased or redundant calibration variants with absolute timing information. Finally, in §6 we described crossfrequency calibration approaches with terms weighted by the frequency-dependent thermal variance σ 2 T (f ). A novel calibration approach could instead apply an ηdependent weighting scheme, where η is the Fourier dual of frequency. One could then downweight delay modes that lie in the 'foreground wedge,' calibrating preferentially to modes that are uncontaminated by sky model errors. 21 cm cosmology is pushing the limits of precision radio interferometry. Progress in the field will require excellent systematic suppression, including the mitigation of calibration errors. Novel calibration approaches can potentially improve measurement sensitivity. These approaches expand upon existing sky-based and redundant calibration techniques, relaxing the assumptions inherent in those approaches and incorporating more physical instrument models.

A. WORKING WITH SINGULAR COVARIANCE MATRICES
In this section we explain how to use Singular Value Decomposition (SVD) to resolve the problem of non-invertable covariance matrices. This approach is fully general and can be applied to any singular covariance matrix, not just those corresponding to traditional redundant calibration as described in §3.2. It remaps an initial calibration parameter basis u init , where the associated covariance matrix C M,init = cov[u init u † init ] is singular, to a new basis u with an invertable covariance matrix C M = cov[uu † ].
All covariance matrices have eigenvalues λ j ≥ 0 and eigenvectors that span the space. We can therefore rotate u init into the basis of eigenvectors of C M,init . We define a new vector u diag as u init expressed in the orthogonal eigenbasis of C M,init : and e j are the orthogonal eigenvectors of C M,init . Here N init is the length of u init . The associated covariance matrix is C M,diag = cov[u diag u † diag ] = EC M,init E −1 = diag( λ 1 λ 2 · · · λ Ninit ), where λ j are the eigenvalues of C M,init . If C M,init is singular than it will have λ j = 0 for some j. This has physical meaning. If λ j = 0 then var[u diag,j ] = 0. It follows that those elements of u diag must take their expectation values: u diag,j = u diag,j . We assume that u init = m init , so likewise u diag = m diag where m diag = Em init .
Therefore u diag,j = m diag,j for all j where λ j = 0. Those calibration parameters simply take on their modeled values; they are no longer fit during the calibration process.
We can now restrict the basis of the calibration parameters to only those that can vary in calibration. The number of independent calibration parameters is N red , which is equal to the rank of C M,init . The new basis of calibration parameters is given by u = Tu init .
Here T is where we use only e j for λ j = 0. C M = cov[uu † ] is now diagonal and invertable.
To develop physical intuition for this remapping we can consider the case of a perfectly redundant array as described in §3.2, where u init correspond to the full set of visibilities. Here cov[u init,j , u init,k ] = var[u init,j ] if baselines j and k belong to the same redundant set and cov[u init,j , u init,k ] = 0 otherwise. In this case, eigenvectors e j for λ j = 0 represent the average visibilities from redundant baseline sets. Eigenvectors e j for λ j = 0 represent visibility differences within redundant baseline sets.
T is a rectangular matrix with linearly independent rows. This means that a right pseudoinverse T p exists such that TT p = 1. However, it does not follow that T p T = 1, so inverting Equation A5 is not straightforward. Rather, we get the somewhat unwieldy expression u init = T p u + j for all λj =0 m diag,j e j |e j | 2 .
This expression simplifies when we assume that m diag,j = 0 when λ j = 0 -in other words, that the model values m init are orthogonal to the null space of C M,init . We then get that This is a physically-motivated assumption. For the case of a perfectly redundant array, requiring that m init is orthogonal to the null space means requiring that model visibilities from different baselines within redundant baseline sets are equal. Components of m init in the null space indicate a disagreement between the instrument models used to produce m init and C M,init .