New Constraints on Lorentz Invariance Violation from Combined Linear and Circular Optical Polarimetry of Extragalactic Sources

Expanding on our prior efforts to search for Lorentz invariance violation (LIV) using the linear optical polarimetry of extragalactic objects, we propose a new method that combines linear and circular polarization measurements. While existing work has focused on the tendency of LIV to reduce the linear polarization degree, this new method additionally takes into account the coupling between photon helicities induced by some models. This coupling can generate circular polarization as light propagates, even if there is no circular polarization at the source. Combining significant detections of linear polarization of light from extragalactic objects with the absence of the detection of circular polarization in most measurements results in significantly tighter constraints regarding LIV. The analysis was carried out in the framework of the Standard-Model Extension (SME), an effective field theory framework to describe the low-energy effects of an underlying fundamental quantum gravity theory. We evaluate the performance of our method by deriving constraints on the mass dimension $d=4$ CPT-even SME coefficients from a small set of archival circular and linear optical polarimetry constraints and compare them to similar constraints derived in previous works with far larger sample sizes and based on linear polarimetry only. The new method yielded constraints that are an order of magnitude tighter even for our modest sample size of 21 objects. Based on the demonstrated gain in constraining power from scarce circular data, we advocate for the need for future extragalactic circular polarization surveys.


Introduction
Einstein's theory of general relativity provides an excellent classical model of gravitation, and the Standard Model of particle physics is a well-established quantum theoretical model of particles and all forces except gravity. Together, they provide a well-tested description of nature at experimentally attainable energies. However, at the Planck scale (E P ≈ 1.22 × 10 19 GeV), a quantum-consistent theory of gravity is required. A lack of direct experimental input to guide the development of the theory poses a significant challenge. Additionally, the failure of the Large Hadron Collider to detect evidence regarding physics beyond the Standard Model, including supersymmetry, presents a challenge to many candidate theories [1] (Indications of beyond-the-standard-model physics at the 3.1σ level in B quark decays were presented by the LHCb collaboration after completion of this work [2]).
Several theories that attempt to unify gravity and the Standard Model at the Planck scale suggest that there may be deviations from Lorentz invariance at this energy scale [3][4][5][6][7][8][9][10]. This motivates detailed tests of Lorentz symmetry despite the fact that such deviations are expected to be highly suppressed at energies E E P [11,12]. Tests of this kind are routinely carried out by high-energy physics experiments [13][14][15][16][17]; however, reaching progressively higher energies will eventually become unfeasible. Astrophysical tests in the photon sector have proven to be particularly powerful because tiny deviations from the speed of light as a result of Lorentz invariance violations (LIV) accumulate when photons propagate over very large distances resulting in potentially observable effects.
If Lorentz symmetry is broken, the phase velocity of light in vacuum may depend on the photon energy, polarization, and direction of propagation. The Standard-Model Extension (SME) is an effective field theory framework that extends the Standard Model of particle physics by introducing new, Lorentz, and CPT violating terms in the Lagrangian, while conserving the charge, energy, and momentum [12,[18][19][20][21]. Within this framework, group theory considerations allow a classification of potential quantum gravity models in three broad classes with respect to their Lorentz violating effects: birefringent and non-birefringent CPT-even models as well as CPT-odd models, all of which result in a birefringent photon dispersion [12]. Non-birefringent models result in a dispersion relation that may depend on the photon energy and propagation direction but will not exhibit any helicity dependence. The strongest constraints on models of this kind result from astrophysical time-of-flight measurements of gamma-rays emitted by transient events or variable sources [22][23][24][25][26][27][28][29][30][31].
Birefringent CPT-even and CPT-odd models can be constrained more strongly by polarization measurements because the measurement is essentially that of a phase difference between the two polarization modes rather than photon arrival times [12,32]. Most astrophysical radiation processes result in very low circular polarization; however, linear polarization can be significant. Birefringence, then, results in a rotation of the linear polarization direction as photons propagate. If the strength of this effect depends on energy, an energy dependence of the polarization direction will be observed, even if the polarization at the source does not depend on energy.
When measuring photon polarization over a broad bandwidth, this rotation results in an effective reduction of the observable polarization fraction. Hence, any observation of linear polarization of light from a distant object can be used to constrain the magnitude of birefringence due to LIV [12]. Some of the strongest constraints of birefringent LIV models come from X-ray polarization measurements [33][34][35][36][37]. However, observations of a single source cannot be used to constrain anisotropic models. In CPT-even models, the LIV terms result in a coupling between the two helicity states, which means neither linear nor circular polarization of light is preserved during propagation. The helicity of ±2 of the coupling necessarily results in anisotropy.
Within the SME, models can be classified based on their low-energy behavior described above, as well as an expansion in terms of the mass dimension d of the corresponding operators, and an expansion in spherical harmonics describing anisotropic effects. The terms of the expansion are then characterized by a set of coefficients, which can be constrained by experiment [38]. In the photon sector, terms of odd mass dimension d represent CPT-odd models and even-d terms represent CPT-even models. Odd-d models are characterized by a set of complex coefficients with (d − 1) 2 real components. Non-birefringent even-d models are characterized by (d − 1) 2 real components, and birefrigent even-d models have 2(d − 1) 2 − 8 coefficients with the same number of real components [12].
We developed a method to combine polarization measurements from many objects in order to individually constrain the coefficients of a given mass dimension. By applying this method to a large number of optical polarization measurements, we obtained the strongest constraints on individual coefficients to date [32,[38][39][40]. In essence, we sample the SME coefficient space and, for each set of coefficients, calculate the likelihood to make all given polarization measurements with the published uncertainties of the linear polarization values. Here, we improve on our approach for CPT-even models characterized by even-d coefficients by incorporating circular polarization measurements of active galactic nuclei in the optical band.
As we will show in Section 3, in these models, the coupling between left-handed and righthanded circular polarization will dominate the reduction of observable linear polarization in most cases. Hence, circular polarization measurements provide an important additional constraint resulting in significantly tighter constraints than with linear polarization alone. In this paper, we extend our method to include circular polarization data. We then apply the method to a set of 21 linear and circular polarization measurements of quasars in order to derive new constraints on the 10 birefringent SME coefficients of mass dimension d = 4. Our new constraints are an order of magnitude stronger than the existing constraints on the individual d = 4 photon-sector SME coefficients.
The remainder of the paper is structured as follows. In Section 2, we briefly describe the photon-sector Lagrangian and photon dispersion relation in the minimal SME, and then derive the underlying equations of our method in Section 3. Section 4 summarizes the assumptions regarding linear and circular polarization at the source that we make in our analysis. The dataset is described in Section 5. We explain the Markov-Chain Monte Carlo method used to sample the SME coefficient space and give the results of our analysis in Section 6. We conclude with a summary in Section 7.

Photon Sector SME
For photons in a vacuum, SME operates by adding two extra terms to the standard Lagrangian of the electromagnetic field, such that the total Lagrangian reads [12,41]: where the added terms are placed in brackets, F αβ is the field tensor, A α is the 4-potential of the field (F αβ = ∂ α A β − ∂ β A α ), ε γδαβ is the Levi-Civita symbol, (k AF ) γ and (k F ) γδαβ are the Lorentz invariance-violating operators, and the sets of coefficients k 4) quantify the effect of the SME. The coefficients are grouped by the mass dimension of the corresponding term in the Lagrangian, d. All coefficients must vanish identically if the standard electromagnetic Lagrangian holds perfectly and no LIV effects are present in the universe. The added terms allow for all possible violations of Lorentz invariance in rotations and boosts of the electromagnetic field, while maintaining Lorentz invariance for the inertial frame of the observer and, thereby, ensuring that physics is independent of the chosen system of coordinates [20]. Specifically, non-zero components of k (d) F γδαβλ 1 ...λ (d−4) give rise to CPTeven terms in the Lagrangian that preserve the CPT symmetry, while non-zero components of 3) result in CPT-odd terms that violate both CPT and Lorentz invariance.
In this work, we follow [20] and further restrict our attention to the so-called minimal SME that only contains terms of renormalizable mass dimensions, i.e., d ≤ 4. This simplification is motivated by the scarcity of the available circular optical polarimetry of extragalactic sources, which necessitates a theory with the smallest number of free parameters in order to derive reliable constraints. We, however, emphasize that the method presented here is universal and can be easily extended to higher mass dimensions (d > 4) if provided with a sufficient amount of experimental data from future polarimetric surveys and necessary computational resources.
Under minimal SME, only the d = 3 term in Equation (2) and d = 4 term in Equation (3) remain, resulting in the following Lagrangian: where k AF has the units of mass and four independent components, and k F is dimensionless and has 4 4 = 256 components, of which only 19 are independent due to the required symmetries: (k (4) [20]. The equations of motion associated with the Lagrangian in Equation (4) are the modified Maxwell equations with SME-induced LIV. Two leading order plain wave solutions (eigenmodes) exist with orthogonal electric field vectors and a phase shift. Following the conventions in [12,41], we write out the modified dispersion relations for the two solutions as follows: where E and p are the energy and momentum of the photon (setting c =h = 1) and ς α are functions of the direction of arrival at the observer and the wavelength of the photon with the functional forms determined by the components of k     AF . The observed electromagnetic wave is a superposition of both eigenmodes.
A non-zero value of ς 0 modifies the speed of propagation (E/p) for both eigenmodes equally and independently of the photon wavelength for renormalizable mass dimensions. Since ς 0 has no effect on polarization, no constraints on its value can be derived in this study. However, ς 0 gains wavelength dependency at d > 4, in which case its effect may be detected through vacuum dispersion (the dependence of the speed of light on the wavelength) [22,41,42]. Non-zero values of ς 1 , ς 2 , and ς 3 result in a relative difference in the propagation speed between the two eigenmodes and will, therefore, alter the polarization state of the observed wave in-flight (vacuum birefringence).

Monochromatic Observations
It is convenient to describe electromagnetic polarization in the Stokes formalism [43] where each polarization state is uniquely identified by a Stokes vector, s s s = (q, u, v), where q, u, and v are the intensity normalized Stokes Q, U, and V parameters. In this formalism, v sets the total fraction of circular polarization (between −1 and +1 with negative and positive fractions corresponding to clockwise and anticlockwise modes respectively, looking outward such that anticlockwise is due East) and q and u encode the linear polarization fraction, p, and angle, ψ, according to the equations below: where atan2(y, x) is the inverse tangent of y/x in the appropriate quadrant (atan2(y, x) = arg(x + iy) for real x, y). Following the International Astronomical Union (IAU) convention ( [44], comm'n 40.8), the frame of reference in the Stokes space is defined such that ψ is measured from North due East in the International Celestial Reference System (ICRS) described in [45,46]. A purely linear polarization state requires v = 0 and p = 0. The eigenmodes in Equation (5) are then given by two antiparallel Stokes vectors (ς 1 , ς 2 , ς 3 ) and (−ς 1 , −ς 2 , −ς 3 ) [47]. Both vectors represent the so-called birefringence axis, which we denote with ς ς ς and, for convenience, take to be positive: The polarization state of a photon (composed of both eigenmodes) propagating in the SME universe, s s s will then precess around the birefringence axis in uniform circular motion with the revolution period of π/(ω|ς|), corresponding to the time taken by the phase shift between the two eigenmodes in Equation (5) to evolve through 2π. The precession direction is against the right hand rule with respect to the birefringence axis. The process is schematically depicted in Figure 1. In the figure, the intensity-normalized Stokes parameters are presented as Cartesian dimensions in 3D space (so-called Stokes space). The locus of all purely linear polarization states is represented by a horizontal plane passing through the origin. Mathematically, the equation of motion is given by: The photon is emitted at the source in an initial polarization state, whose location in the Stokes space (shown here) is indicated with vector s s s z . In-flight, the state will precess around the birefringence axis ς ς ς through angle 2Φ until, eventually, the photon arrives at the telescope in the state s s s 0 . The direction of the birefringence axis and the rate of precession are determined by the particular SME configuration. The blue quadrilateral represents the plane of linear polarization where Stokes V is 0. Individual components of ς ς ς = (ς 1 , ς 2 , ς 3 ) and the angles referenced in text (ξ, ζ) are labeled as well.
In Equation (9), the factor of 2 arises from the fact that the propagation speed between the two eigenmodes differs by 2|ς| as per Equation (5). We emphasize that the birefringence axis is defined in Stokes rather than physical space, and, as such, its existence does not immediately imply a preferred direction in the universe. Instead, the anisotropy is generated by the directional dependence of the birefringence axis discussed in Section 3.2.
The Stokes vector represents the polarization state of a photon arriving from a specific direction in the sky and is, therefore, a function on the celestial sphere, requiring a local reference direction of ψ = 0 (ICRS North according to the IAU convention). v is independent of the choice of reference and q ∓ iu rotate as spin ±2 objects. In this context, a function f of celestial coordinates (such as q, u, v, or q ± iu) is said to have spin s if it transforms as f → f exp(isδψ) when the local North is rotated by δψ radians due East. In this case, Taking advantage of those symmetries, we re-express ς ς ς and s s s in a so-called spin-weighted basis, such that each dimension of the vector space has a defined spin. In particular, we follow the convention used in [12] and adopt the new basis as s s s = (q − iu, v, q + iu) and where the spins of vector components are +2, 0, and −2 respectively. The transformation into the new basis may be written as: This, in turn, re-organizes the equation of motion, Equation (9), into: From this point onward, the spin-weighted basis in Equation (11) will be assumed in all vectors and matrices. Denoting the total precession angle from emission to observation of the photon polarization state in Stokes space around the birefringence axis with 2Φ, the solution to Equation (11) is given by the precession matrix obtained using Rodrigues' rotation formula ( [48], p. 209) and transformed into the spin-weighted space withŜ: whereς ς ς is the normalized ς ς ς and the following definitions are introduced: The angles ξ and ζ have geometric significance and are labeled in Figure 1. In Equation (  12), s s s 0 = (q 0 − iu 0 , v 0 , q 0 + iu 0 ) is the observed Stokes vector of the photon on Earth, and s s s z = (q z − iu z , v z , q z + iu z ) is the initial Stokes vector at redshift z.
Considering the precession period of π/(ω|ς|) for the polarization vector, the equation of motion for Φ reads: where T is the total time of flight for the photon, z is the redshift of the source, and ω 0 is the observed angular frequency of the photon at z = 0. We used dt = −dz/[(1 + z)H z ] and ω(z) = ω 0 (1 + z). The Hubble expansion parameter, H z , is given by In Equation (17), H 0 is the Hubble constant and Ω r , Ω m , Ω k , and Ω Λ are the present day fractional contributions of radiation, matter, curvature, and dark energy to the energy budget of the universe. In this study, we adopt H 0 = 67.66 km s −1 Mpc −1 = 1.4433 × 10 −33 eV, Ω m = 0.3111, Ω r = 9.182 × 10 −5 , Ω Λ = 0.6889, and Ω k = 0, following [39,49] (see [39] in particular for a discussion of the accuracy of those values in this context and the Hubble tension).
While, in principle, both CPT-even (d = 4) and CPT-odd (d = 3) contributions to the Lagrangian may exist simultaneously, it is generally preferable to begin the search for Lorentz invariance violation with the smallest number of free parameters. We will follow [12,22,32,39,40,50] and assume that the effects at a specific d dominate. In particular, we choose to focus on the d = 4 case as an even mass dimension is required for helicity coupling (for odd d, the birefringence axis is parallel to the Stokes V axis and, therefore, circular polarization of the photon is conserved). Furthermore, since ς 3 ∝ ω −1 ∝ ω −1 0 (see Equation (2)), Φ is uniquely wavelength-independent in the d = 3 case, which motivates alternative approaches in [51][52][53][54][55]  AF vanish and, therefore, so does ς 3 . Geometrically, the d = 4 case has the birefringence axis, ς ς ς, constrained to the plane of linear polarization. For, ς 3 = 0, ζ (Equation 14) evaluates to π/2, allowing the precession formula in Equation (12) to be simplified as follows:

Directional Dependence
In Equation (18), the dependence of Φ on the direction of arrival is determined by the particular SME configuration, which we seek to constrain by comparison with observations. Following [12,39,41], we limit the dependence to 10 linearly independent degrees of freedom (corresponding to the 10 independent components of k (4) F setting ς ± ) by expanding ς ± in terms of l = 2, s = ±2 spin-weighted spherical harmonics, s Y l,m , where s, l, and m are the spin, azimuthal, and magnetic numbers, respectively, (−2 ≤ m ≤ 2): Here,n n n is a unit vector pointing in the direction of photon arrival (opposite to the photon momentum) and k (E,B)l,m are complex coefficients of the expansion that obey the following symmetries: such that there are 10 real independent parameters: k (E)2,0 , k (B)2,0 , Re k (E)2,1 , Im k (E)2,1 , Re k (E)2,2 , Im k (E)2,2 , Re k (B)2,1 , Im k (B)2,1 , Re k (B)2,2 , and Im k (B)2,2 .
Since none of the terms in Equation (19) have spherical symmetry, the d = 4 case considered in this study is naturally anisotropic. However, isotropic terms do exist in odd d SME configurations [12,22,32,39].
For astronomical tests,n n n may be written in some system of celestial coordinates. In this work, we usen n n = (α, δ), where α and δ are the ICRS [45,46] right ascension and declination.
Substituting Equation (19) in Equation (16): Φ is directly proportional to the comoving distance to the photon source ( z 0 dz /H z ), which is a unique feature of d = 4 SME.

Broadband Observations
Astronomical observations of the Stokes parameters q 0 , u 0 , and v 0 are taken over some finite band defined by a dimensionless detection efficiency τ(ω 0 ), rather than a single monochro-matic frequency. The observed effective Stokes parameter, Q 0 , is then determined by the weighted average of the observed Stokes parameter q 0 : where we define the operator T as the efficiency-weighted average. The other two Stokes parameters are then The observed broadband linear polarization fraction (Π 0 ) and angle (Ψ 0 ) may then be defined equivalently to Equations (6) and (7): For convenience, we introduce a primed system of Stokes coordinates, rotated with respect to the ICRS coordinates anticlockwise through ξ: As a result of this rotation, the polarization angle decreases by ξ/2, while both the linear and circular polarization fractions remain unchanged. In primed coordinates, the observed broadband Stokes parameters can be expressed in simple forms by assuming that the source polarization spectrum only slowly varies with ω across the optical range, implying that s s s z (ω z ) ≈ s s s z . This assumption is justifiable by the low redshift (i.e., SME-free) spectropolarimetry [32] that shows a weak dependence of polarization on wavelength for most galaxies due to the optical range being comparatively narrow. Q 0 , U 0 , and V 0 are then related to the Stokes parameters at the source as: The functional forms of T [cos(2Φ)] and T [sin(2Φ)] are entirely determined by the instrument and may be pre-tabulated for a range of Φ values for fast lookup. Equations (26)- (28) demonstrate that the total polarization fraction ( q 2 + u 2 + v 2 ) is conserved in the monochromatic limit (T [x] → x); however, it will, in general, decrease with redshift for broadband observations due to the polarization washout induced by the frequency dependence of d = 4 SME effects (i.e., |T [x]| < |x|).
This effect is illustrated in Figure 2, displaying the value of T [sin(2Φ)] as a function of Φ for the cases of monochromatic, filtered broadband and unfiltered broadband observations. While the amplitude of the operand is conserved in the monochromatic case, it decays in both broadband cases with the unfiltered decay being more rapid due to a wider range of wavelengths observed.   (27) and (28)) as a function of Φ (see Figure 1). The three cases shown correspond to a monochromatic observation (black, T [sin(2Φ)] → sin(2Φ)), a broadband observation with an unfiltered Ga-As photomultiplier tube (red), and a broadband observation in the V-band of the EFOSC2 [56] instrument (green). This calculation assumes observations at the zenith. See Section 5 for more details on instruments, bands, and atmospheric effects.

Likelihood Model
To differentiate SME predictions for the observed broadband polarization parameters (Q 0 , U 0 , V 0 , Ψ 0 , and Π 0 ) derived above from the actual experimental measurements, we introduce subscript m for the latter, such that, in the case of the chosen SME configuration being a perfect description of LIV, we must have Q 0 = Q m , U 0 = U m etc. in the absence of other effects.
In the limit of weak SME motivated by the strict constraints established in previous works [32,39,40], we expect the change in Stokes parameters due to SME to be small, such that the overall ratio of linear to circular polarization (Π/V) for extragalactic photons is only mildly perturbed in-flight. Since both experimental and theoretical considerations (e.g., [57][58][59][60]) suggest that linearly polarized emissions dominate over circularly polarized emissions for nearly all realistic extragalactic sources, we further assume that |v z | p z and |V 0 | Π 0 . Under this assumption, Equations (27) and (28) suggest that the overall effect of the SME is to suppress linear polarization and enhance circular polarization in-flight, as demonstrated in Figure 3b where linear and circular polarization fractions are plotted on the same set of axes as functions of redshift. Therefore, SME coefficients are constrained by measurements of lower |V m | and higher Π m . In the extreme case, an error-free measurement of V m = 0 or Π m = 1 would entirely rule out all SME configurations considered in this study.
Assuming that the source polarization is somehow known (see Section 4), Equations (26)-(28) can be used to predict the expected observed polarization for any d = 4 SME configuration and for any instrument (encoded in the functional form of T ). The prediction must then be compared to experimental results, which may be accomplished by evaluating the likelihood of compatibility of the observed data with theory. For broadband polarimetry, three types of measurements are available: • linear polarization fraction, Π m , • polarization angle, Ψ m , and • circular polarization fraction, V m .
In our method, we advocate deriving constraints on SME parameters from the linear and circular polarization fraction measurements, Π m and V m . The polarization angle, Ψ m , is still essential in the method (Section 4), but will not be used in the likelihood model directly.
Since the measured circular polarization, V m , may take both positive and negative values, we assume a Gaussian parent distribution with the mean of V m and the standard deviation of σ V . The likelihood of compatibility of such a measurement with the SME prediction V 0 may be written as follows: The integration bounds over the Gaussian distribution in the equation above are chosen such that a measurement (V m ) is considered compatible when its absolute value (|V m |) is as large or larger than the SME prediction (|V 0 |). In other words, P circ is the probability of The measured linear polarization degree, Π m , may display strongly non-Gaussian behavior due to being an intrinsically positive quantity. Following the derivation in Appendix A, we assume Π m to be drawn from the Rice distribution (P(Π|Π), Equation (A5)) withΠ = Π m and Var[Π] = σ 2 Π , where σ Π is the uncertainty of the measurement and Var[Π] is given in Equation (A7). The statistical distribution of the linear polarization fraction and the bias in linear polarization measurements are discussed in Appendix A.
Since the overall SME effect is to decrease the linear polarization, a particular measurement (Π m ) is considered compatible with a particular SME configuration if it is smaller than the SME prediction (Π m < Π 0 ). Therefore, the probability of compatibility with the SME prediction Π 0 may be written as follows: P lin = Π 0 0 P(Π|Π m )dp (30) where P(Π|Π m ) is the Rice distribution from Equation (A5). The total compatibility of a given dataset of both circular and linear polarimetry is obtained by multiplying all individual likelihoods together: where the product is taken over all measured sources. The constraints on the 10 independent SME parameters are derived by maximizing the likelihood with respect to them and extracting the standard errors, as detailed in Section 6.

Source Parameters
The likelihood model introduced in the previous section requires some assumption of the photon polarization state at the source, s s s z . In this section, we introduce our assumptions for the polarization angle and linear and circular polarization fractions at the source: ψ z , p z , and v z . In general, we aim to choose the most conservative values, i.e., values for which fewest SME configurations are ruled out or, alternatively, values that maximize the compatibility likelihood P. By maximizing P, we seek the weakest possible constraints on LIV supported by the data. We emphasize that weaker constraints do not imply that a future positive detection is more likely. Instead, they merely define the region of parameter space that we can claim to have thoroughly explored with high confidence.

Circular Polarization
For all observed astronomical sources, we assume the circular polarization at the source to be zero (v z = 0). This assumption is justified by both experimental and theoretical studies in literature [57][58][59][60] that found negligible or very small circular polarization for the majority of extragalactic sources. For sources with small but non-negligible circular polarization (0 < |v z | p z ), the assumption of v z = 0 remains conservative, since non-zero v z results in an increase in the magnitude of observed circular polarization (Equation 28) and, therefore, weaker constraints on SME coefficients (larger P circ ).

Linear Polarization
In our previous works [32,39,40], large values of p z were adopted as most conservative since P lin increases with p z in the limit of weak SME (small Φ) and the dominance of linear polarization over circular (|v z | p z ). Specifically, the most conservative assumption for the value of p z was taken as the largest possible source polarization fraction. For example, the authors in [32,40] adopted p z = 1.0 for all sources ignoring astrophysical considerations, while [39] adopted p z = 0.7 based on the literature polarimetry of low redshift sources, resulting in somewhat stricter constraints on SME parameters.
However, previous work disregarded the fact that P circ , in general, decreases with p z since larger linear polarization at the source results in the faster conversion of linear to circular polarization in-flight and, hence, larger observed circular polarization V 0 . Therefore, the assumption of large p z may not be most conservative in certain instances when circular polarimetry is available. An example of this effect is demonstrated in Figure 3a.
In the figure, the values of P lin , P circ , and the total compatibility P are shown for a test observation and a test SME configuration as a function of the assumed source polarization p z . P lin increases with p z monotonically with an asymptote at P lin = 1, suggesting higher values of p z as most conservative in the absence of other considerations. On the other hand, P circ decreases with p z monotonically. Therefore, a more conservative assumption for p z is one that maximizes the total compatibility, which, in this example, falls around p z = 0.55.
In this study, we determined the most conservative value of p z by numerically maximizing the total likelihood of compatibility, P lin P circ , for each available polarization measurement.

Polarization Angle
The polarization angle at the source, ψ z , may be estimated by seeking a value that reproduces the observed polarization angle, Ψ 0 = Ψ m , under the most conservative assumptions for p z and v z as described above. Combining Equations (26)-(28) and solving for ψ z : where Q m = Π m cos(2Ψ m ), U m = Π m sin(2Ψ m ), and Ψ m = Ψ m − ξ/2. As in [39,40], the uncertainty in Ψ m is not required by the method. The astrophysical source is assumed to be located at the Vernal equinox (α = 0, δ = 0) and z = 2.0. For demonstration purposes, the adopted test measurement is Π m = 0.5 ± 0.1, V m = 0.00 ± 0.01 and Ψ m = 0.0. The adopted SME configuration has all ten parameters considered in this study set to 10 −34 (the order of magnitude of the upper limit derived in [39,40]). The test measurement is assumed to have been taken through the EFOSC2 [56] V filter in zenith (see Section 5). The linear compatibility increases while the circular compatibility decreases with p z . The most conservative assumption for the value of p z is the one maximizing the total compatibility, which, in this case, occurs around p z = 0.55. (b) The effect of SME on the linear and circular polarization fractions, shown here as Π 0 and V 0 as functions of the redshift of the source. The predictions shown here were derived for a test source at the Vernal equinox (α = 0, δ = 0) observed through the EFOSC2 [56] V filter in zenith. The initial polarization state of the photons was assumed to be v z = 0, p z = 1, and ψ z = −30 • . The adopted SME configuration is the same as in (a).

Sample Dataset
We derive constraints on the 10 real SME parameters by maximizing the total likelihood of compatibility (Equation 31) for a set of linear and circular polarimetric measurements from literature. The dataset used in this study includes 21 quasars and was adopted from the compilation of archival linear polarimetry and original circular polarimetry in [57]. All data used in this analysis are reproduced in Table 1.  .) are provided for the linear polarization measurements. All circular polarization measurements were taken through the Bessel V filter and were reported in [57]. The detect efficiency profiles for the Bessel V filter, Ga-As photomultiplier, and Na-K-Cs-Sb (S20) photomultiplier are shown in Figure 4. Circular polarization measurements in [57] were obtained through the Bessel V filter (ESO #641) using the EFOSC2 [56] instrument mounted on the ESO 3.6 m telescope at La Silla Observatory.
The measurements from Impey+1990 [66] were obtained with an unfiltered Ga-As photomultiplier tube in the MINIPOL polarimeter mounted on the Irénée du Pont 100 inch telescope at the Las Campanas Observatory. The measurements from Berriman+1990 [73] were obtained with an unfiltered Ga-As photomultiplier in the Two-Holer polarimeter/photometer mounted on the Bok 2.3 m telescope at the Steward Observatory and the 1.5 m telescope at the Mount Lemmon Observing Facility. The measurements from Schmidt+1999 [72] were obtained unfiltered using polarimeters mounted on the Bok 2.3 m telescope at the Steward Observatory. The publication does not specify the detector; however, it is mentioned that the waveband and sensitivity were very similar to the Breger polarimeter [79] at the McDonald Observatory.
Since [79] employ the response curve of a Ga-As photomultiplier tube in their analysis, we assume that the Steward Observatory polarimeters used some type of Ga-As tubes as well. The measurements from Visvanathan+1998 [78] were obtained with an unfiltered Na-K-Cs-Sb (S20) photomultiplier in an automated polarimeter [80] mounted on the Anglo-Australian Telescope at the Siding Spring Observatory. The measurements from Hutsemékers+1998 [62], Lamy+2000 [68], and Sluse+2005 [64] were obtained through the Bessel V filter using the EFOSC1/EFOSC2 [56] instrument at the La Silla Observatory, which we assume to have a sufficiently similar transmission profile to the Bessel V filter (ESO #641) used by [57] for circular polarimetry.

Detection Efficiency Profiles
Our method requires efficiency profiles, τ(ω 0 ), for each measurement in the dataset. For filtered observations, the transmission profile of the filter is typically self-sufficient with the atmospheric transmission, detector response, mirror reflectivity, etc. only contributing minor higher order corrections. For unfiltered observations, the response curve of the detector becomes the most determining factor, with the atmospheric transmission making a significant contribution at the blue end of the visible spectrum where the ozone layer predominantly sets the short wavelength cut-off of the instrument.
For all filtered observations through the Bessel V band, we adopted the ESO #641 transmission profile available online (https://www.eso.org/sci/facilities/lasilla/instruments/ efosc/inst/Efosc2Filters.html, accessed "13 May 2021"). For unfiltered observations, we took the generic response curves of the corresponding photomultiplier tubes (Ga-As and S20) from ( [81], p.101). All profiles were multiplied by the atmospheric transmission, based on the atmospheric radiation model in [82] calculated for the geographic altitude of 2400 m and weather conditions of Cerro Paranal. The model employs the radiative transfer routine in [83] and HITRAN opacity database [84]. The altitude and conditions adopted by the model are generally similar to those at the majority of aforementioned facilities contributing data to this study.
Since the atmospheric transmission strongly depends on the target elevation above the horizon, we performed all calculations for the cases of both optimal (airmass Z = 1.0) and poor (airmass Z = 3.0) conditions. In this analysis, we define the airmass, Z, as the integrated atmospheric density along the line of sight to the source normalized to Z = 1 in zenith. The final efficiency profiles used in both cases are provided in Figure 4. All profiles are normalized to τ = 1 at the peak efficiency for convenience, since our method is not sensitive to multiplicative pre-factors, as evident from Equation (22).

SME Constraints
Following [39,40], we explored the ten-dimensional likelihood distribution given by Equation (31) using the Metropolis-Hastings Markov-Chain Monte Carlo (MCMC) approach implemented in Python package emcee [85]. Exploration was carried out by so-called walkers that spawn at some location in the likelihood space (i.e., some combination of SME parameters) and proceed in random Monte Carlo steps. In each step, offsets for the positions of the walkers were drawn from a suitable proposal distribution in all ten dimensions, and the likelihood of compatibility at the new location was estimated.
The ratio of the new likelihood to the old was then compared to a uniform random number between 0 and 1, and the step was accepted if the former exceeded the latter. The nonzero probability of stepping into a lower likelihood was implemented to prevent walkers from becoming "stuck" in local minima if any happened to be present in the likelihood space. The probability distribution for each of the SME parameters was then extracted by assuming the probability of each value along a given dimension to be proportional to the number of steps that the walkers spent in its immediate proximity.
In this study, we followed [39,40] and used a Gaussian proposal distribution with the same standard deviation in all dimensions. The standard deviation was chosen by first assuming all SME parameters to be equal and searching a value where the likelihood was evaluated to the midpoint between the extreme cases of infinitely strong (k (E,B)l,m → ∞) and infinitely weak (k (E,B)l,m → 0) SME. The result, ∼ 2 × 10 −36 , was expected to be of comparable magnitude to the average width of the likelihood distribution across all dimensions. The initial positions for 20 MCMC walkers were drawn from the proposal distribution as well.
The run was terminated when each walker completed 500 steps (accepted or rejected) for a total of 10 4 steps across all walkers. At the end of the chain, the step acceptance ratios were found to be 0.21 and 0.23 for the Z = 1 and Z = 3 cases. The extracted probability distributions for all 10 SME parameters are given in Figure 5. The constraints on the individual parameters (upper and lower bounds) were taken as the 5th and 95th percentiles of the corresponding distributions. The numeric values are summarized in Table 2. In the case of differing results for the Z = 1 and Z = 3 cases, the most conservative case was taken (i.e., the minimum value for the lower bound and maximum value for the upper bound). Table 2. The derived constraints on the values of the 10 real d = 4 SME parameters. The upper and lower bounds were taken as the 5th and 95th percentiles of the distributions shown in Figure 5. For each parameter, the most conservative of the Z = 1 and Z = 3 cases is shown.

Conclusions
In this work, we used optical circular and linear polarimetry of high redshift extragalactic sources to derive new constraints on the Lorentz invariance violation in the framework of the Standard-Model Extension at mass dimension d = 4. The numerical values of the constraints are expressed as upper and lower bounds on the 10 SME parameters summarized in Table  2. The detailed probability distributions for each parameter were also calculated and are presented in Figure 5. Since the exact conditions of the literature observations employed in this study are unknown, all calculations were performed for two different airmasses (Z = 1 and Z = 3) with no significant difference in the results.
Our results assume that SME effects at mass dimension d = 4 dominate. The focus on this particular mass dimension is justified by the fact that it represents the leading-order terms that can be tested with our method. Furthermore, the reliance of our method on helicity coupling, which only occurs at even d, and the large number of independent SME parameters at d ≥ 6 exceeding the number of lines of sight available in our dataset (21) prevented application of our dataset to higher mass dimensions.
For example, d = 6 analysis requires at least 42 lines of sight, d = 8 requires 90 lines of sight, etc. Nonetheless, the method proposed in this work can be extended to those mass dimensions given a larger number of measurements. A follow-up study of higher mass dimensions would ideally also use higher-energy data due to the stronger suppression of these terms. However, circular polarization measurements of astrophysical X-rays are currently unavailable.
In our analysis, we assumed the dominance of linear polarization over circular at all sources as well as sufficiently weak LIV such that said dominance was maintained throughout the line of sight. Both assumptions are well justified by theoretical and experimental astrophysical considerations [57][58][59][60] as well as existing SME constraints [38][39][40]. Due to the difficulty of determining the initial polarization state of photons at the source, we searched the space of all possible initial polarization states numerically and assumed the one that resulted in the most conservative SME constraints for any given source. Finally, we assumed a weak wavelength dependence of polarization at the source across the optical range justified by the nearly flat polarization spectra of nearby active galactic nuclei as reported in our previous work [32].
Similarly to [38][39][40], we found the probability distributions for individual SME parameters ( Figure 5) were consistent with zero within two sigma. Unlike [38][39][40], our results display more asymmetry around the origin, with the most notable example being Im k (E)2,1 , which appears heavily skewed toward negative values. We note that circular polarimetry is generally expected to be more sensitive to the sign of coefficients than linear polarimetry due to the fact that both clockwise and anticlockwise circular polarization may be induced by the SME. The observed asymmetry in the distribution is likely caused by the particular combination of clockwise/anticlockwise circular polarization measurements along the lines of sight considered in the sample dataset.
The discussion of systematic errors in [39] is mostly applicable to this study as well. In our method, it is assumed that any drift in the photon polarization state between the source and the telescope is induced by the LIV and not by astrophysical processes. We first note that any process that reduces linear or enhances circular polarization would weaken our constraints due to the conservative likelihood model used.
While modelling such effects may improve the derived constraints further, we are at no risk of falsely ruling out viable SME configurations. If unaccounted, said astrophysical processes may eventually establish a lower bound on the parameter constraints derivable using our method. Given the considerable improvement of the constraints derived in this study compared to its predecessors (see below), it is reasonable to assume that such a lower bound has not yet been reached.
On the other hand, astrophysical processes that enhance linear or reduce circular polarization are of much greater concern as they would falsely tighten our constraints on SME effects. Fortunately, very few such processes are known to occur in the optical regime, and most are expected to average out over multiple lines of sight. A prominent example of such an effect includes the polarization and de-polarization of radiation by the interstellar medium [86,87].
Our constraints are directly comparable to those derived in [39,40] (also listed in [38]) and are tighter than both by approximately an order of magnitude due to our use of circular polarimetry, which has not been previously applied in similar studies. We note further that this result was achieved with only 21 unique lines of sight, which is less than the dataset in [40] containing 27 lines of sight and much less than the dataset in [39] containing 1278 lines of sight. Our analysis demonstrates that deriving constraints from circular and linear polarimetry simultaneously is a more efficient method than the methods previously employed in literature as not only does it provide better constraints for fewer sources but it is also free of the somewhat arbitrary assumption of an excessively high initial linear polarization fraction.
We expect that a larger sample size as well as higher quality circular polarimetry may significantly improve the constraints derived in this work. Unfortunately, the scarcity of circular measurements for extragalactic sources, in part due to their characteristically weak signal, imposes strict limits on the maximum improvement one may expect from archival data and calls for new optical polarization surveys. Additionally, the process of estimating the initial linear polarization at the source, p z , through numerical optimization is far more computationally demanding than the methods in [39,40], limiting the number of MCMC steps used in this study (10 4 compared to 0.5 × 10 6 in [39] and 10 7 in [40]) and potentially requiring high performance computing for the adequate processing of a larger dataset in the future.
Author Contributions: R.G. led the analysis of linear and circular polarization data. P.B. collected the dataset of linear and circular polarization data, filter transmission curves, and PMT efficiencies. F.K. co-developed the analysis method. All authors have read and agreed to the published version of the manuscript.
Funding: FK acknowledges NASA support under grant 80NSSC18K0264.

Acknowledgments:
The authors would like to thank Alan Kostelecký, Matthew Mewes, David Mattingly, Henric Krawczynski, Jim Buckley, Floyd Stecker, and Brian Keating for fruitful discussions. We are grateful to the late Andy Friedman without whom this collaboration would never have happened. This work made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Polarization Statistics
In this appendix, we present our treatment of the linear polarization degree distribution and the linear polarization bias. Our argument is an adaptation of similar treatments in the context of X-Ray polarimetry [88] and radio interferometry [89] to the optical regime. Assuming no variability in the source, there exists some "true" linear polarization degree,p, which is related to the "true" intensity-normalized Stokesq andû parameters: p = q 2 +û 2 . (A1) We assume the measured Stokes parameters, q and u, to be normally distributed around their "true" counterparts,q andû, with identical standard deviations of σ. Furthermore, we assume no correlation between q and u. The probability density of observing a particular pair of q and u is then a product of their respective Gaussian distributions: P(q, u|q,û) = 1 2πσ 2 exp − (q −q) 2 2σ 2 exp − (u −û) 2 2σ 2 . (A2) Equivalently, in polar coordinates: where we substituted q = p cos(2ψ), u = p sin(2ψ) and their "true" counterparts. Note the added pre-factor of 2p due to the coordinate transformation, dqdu = 2pdpdψ. Now integrate Equation (A3) over all ψ to obtain the distribution of p independently of ψ: P(p|p) = π 0 P(p, ψ|p,ψ)dpdψ = p σ 2 exp − p 2 +p 2 2σ 2 I 0 pp where I n (x) is the modified Bessel function of the first kind. Rewriting Equation (A4) in terms of the exponentially scaled Bessel function, I 0 (x) = exp(|x|)i 0 (x), alleviates numerical issues at small σ due to diverging I 0 (x): The distribution is plotted in Figure A1. The value of σ in Equations (A4) and (A5) may not be known a priori. To address this issue, we calculate the expected value, E[p], and the variance, Var[p], of p using the derived distribution: 2σ exp −p 2 4σ 2 (p 2 + 2σ 2 )I 0 p 2 4σ 2 +p 2 I 1 p 2 4σ 2 (A6) Note that, in general, E[p] =p due to the polarization degree bias. This effect is demonstrated in Figure A1. For a given polarization measurement (E[p] orp depending on whether the value was debiased or not) and its uncertainty ( Var