Cosmographic constraints on a Gödel-type rotating universe

Rotation is an ubiquitous phenomena in any astronomical setting. Almost all heavenly bodies have some degree of rotation/spin, from the slow rotation of asteroids to the relativistic speeds at the center of galaxies. It is then natural to ponder over the question of whether the universe itself is rotating or not. Could the universe, on the largest scales, be rotating? Or whether the widespread phenomena of rotation is confined to regions where the universe is inhomogeneous? The question of cosmic rotation and its potential role in the structure and evolution of the universe has intrigued physicists for nearly a century. Whittaker [1] and Lemaître [2] contemplated on the nature of a rotating “primeval atom” from which the universe supposedly emerged in relation to the formation of galaxies and large-scale structures.

The study of rotational effects in cosmology started with the work of Lanczos, who modeled the universe as a rigidly rotating dust cylinder with infinite radius [3], and later in the works of George Gamow [4]. However, a more realistic metric with rotation as a solution to Einstein’s field equations was proposed by Gödel in 1949 [5], that is given by the line element,

where ‘$a$’ is a positive number. Gödel’s universe assumes dust as matter with an energy density $\varepsilon$ and a negative cosmological constant ‘$\Lambda$’. In this model, the angular velocity parameter ‘$\omega$’ that quantifies cosmic rotation is given by $\omega^{2}=1/2a^{2}=4\pi G\varepsilon=-\Lambda$. This model provided a theoretical framework for studying rotating cosmologies for many years. Earlier model of Lanczos’ was less physical compared to the Gödel’s model [6].

The metric in Eq. (1) is stationary and incorporates non-trivial global properties of rotating space-times. However, one fundamental problem with this solution was the existence of closed time-like curves (CTC) proved by Gödel himself [7]. Gödel attempted to derive completely causal rotating cosmologies in his last work on rotation and presented a generalization of his original stationary model incorporating expansion, but without proof. First explicit solutions were reported in Refs. [8, 9]. Maitra [8] extended the van Stockum’s solution [6] with rotation as well as shear, and also introduced a straightforward condition for determining the absence of closed time-like curves, ensuring that the model avoids causality violations. Later, multiple solutions of rotating models with or without expansion have been developed, each offering unique insights into the role of rotation in cosmology. Theoretical advancements have given rotating cosmologies some serious considerations, including Einstein-Cartan theory that incorporates spin and torsion. Exact rotating solutions in Einstein-Cartan theory of gravity with spin and torsion were described, for example, in Refs. [10, 11, 12, 13, 14, 15, 16, 17, 18, 19].

The discovery of the cosmic microwave background (CMB) and its high degree of isotropy provided a very strong case for isotropic models. However, deviations from a near perfect isotropy are indicative of possible anisotropies present in the early stages of the evolution of our universe. This motivated a class of homogeneous but anisotropic cosmologies known as Bianchi type cosmological models [20, 21, 22]. Gödel’s solution is one such solution within the class of Bianchi type-III cosmological models. Homogeneous but anisotropic rotating cosmologies were verily analyzed earlier in the literature, in which upper limits on the value of cosmic rotation were obtained from data [23, 24, 25, 26, 27, 28]. But, those studies did not separate the effects of vorticity and shear. In Refs. [29, 30], the effects of pure cosmic rotation were delineated from cosmic shear.

Even with the plethora of astrophysical data sets that were accumulated over the years, no completely convincing results were reported. Birch [31, 32, 33] observed an apparent anisotropy in the alignment between the polarization vectors and the major axes of radio sources. He interpreted this observed anisotropy as potential evidence for cosmic rotation. With the advent of precision era in cosmology, primarily driven by CMB observations, the Cosmological principle which is a foundational assumption of modern cosmology was put to test. This lead to finding many instances of isotropy violation in various astronomical and cosmological data sets (see for example Ref. [34]). Such instances of isotropy violation seen in CMB data have come to be known as CMB anomalies [35, 36, 37, 38, 39]. There were many coincidental alignments among various preferred axes seen in these diverse data sets [40]. Besides isotropy violations, the standard concordance model is also riddled with many tensions [41, 42, 43] from observations that challenge it despite its broad success. Apart from providing possible explanation to various signals of anisotropy seen in data, rotating models may also offer mild corrections to the Hubble flow, potentially providing an alternative explanation for dark energy and a possible resolution to the Hubble tension problem (for a review of Hubble tension and various possible (re)solutions see, for example, Ref. [44, 45, 46, 47]).

A study of spins in the intergalactic medium suggests that longest substantially rotating objects in the Universe are likely filaments of matter connecting neighboring galaxy pairs [48]. Another study reported observational evidence that cosmic filaments spin [49]. It was also reported that there is an asymmetry in the distribution of galaxy spin directions (handness) [50, 51, 52]. Further, some spectroscopic surveys indicate that the direction of rotation of a galaxy is coherent with the average motion of its neighbouring galaxies over a scale of 1 Mpc (megaparsec) [53, 54]. It was argued earlier that primordial turbulence or cosmic vorticity could have played a significant role in structure formation, providing a natural origin for galactic angular momentum [55, 56]. The availability of high-resolution CMB data in the last two decades saw a renewed interest in homogeneous but anisotropic cosmologies, including rotating models [57, 58]. Consequently, more stringent limits on cosmic shear and vorticity were derived using CMB [59, 60, 61].

In Sec. 2, we layout the general properties of a rotating and expanding form of Gödel metric given in Eq. (1) and present the relation between apparent magnitude of a luminous source ‘$m$’ vs redshift ‘$Z$’ in such a model. In Sec. 3, the data used in constraining the Gödel-like model under study and the methodology to derive constraints on its parameters are described. Finally, we present our results in Sec. 4 and summarize our finding in Sec. 5.

In this section we briefly discuss the Gödel-like model that we are going to investigate in this work. First we present an overview of Gödel-type metric with additional details as necessary, followed by a description of the Taylor expanded apparent magnitude vs redshift relation in such a universe. For more details the reader may consult some of the original works and the references therein [29, 62, 63, 64].

The Pantheon+ Type Ia supernova (SNIa) dataset¹¹1https://github.com/PantheonPlusSH0ES/DataRelease contains 1701 light curves corresponding to approximately 1550 unique spectroscopically confirmed SNIa [66]. The remaining light curves in the dataset pertain to either the same supernova (SN) observed by different surveys or to “SN siblings” that are SNIa found within the same host galaxy. This latest data set includes significant enhancements, such as the addition of new SNIa and corrections for peculiar velocity effects and host galaxy bias in the SNIa covariance matrix.

In order to measure the luminosity distance to these supernovae (SNe), we use Cepheid variable stars, a type of pulsating stars whose intrinsic luminosity is related to their period of variability. This relationship, known as the period-luminosity (P-L) relation, was first discovered by Henrietta Swan Leavitt in 1912 [67, 68] while studying stars in the Magellanic Clouds. Cepheids are bright enough to be observed in relatively nearby galaxies. Distances to them within the Milky Way were first measured using parallax (a geometric method). This provides a way to estimate the absolute brightness (or luminosity) of these stars based solely on their observed period. Next, Cepheids in nearby galaxies - those that also host Type Ia supernovae (SNe Ia) - are observed to link the supernova luminosity scale to absolute distances. This step is crucial because once the absolute magnitudes of SNe are calibrated using Cepheid distances, SNe Ia can then be used to measure much greater cosmological distances. This calibration process allows us to determine the luminosity distance to very distant supernovae in the Hubble flow (of the expanding universe). The Pantheon+ dataset contains Cepheid anchors that were part of the Pantheon+SH0ES analysis [69]. 42 SNe Ia from Cepheid-calibrated host galaxies are included in the analysis, along with 277 Hubble-flow supernovae.

SNIa cosmology relies on comparing the theoretical relationship between redshift ($Z^{*}$) and luminosity distance ($\mathcal{D}$) or equivalently the distance modulus ($\mu^{\rm th}$) that is dependent on various cosmological parameters of our model, with the observed distance modulus ($\mu^{\rm obs}$) to derive constraints on them. The theoretical and observed distance modulus can be formally expressed as

where ‘$Z^{*}$’ is the redshift of an SNIa, ${\boldsymbol{\xi}}=\{h_{0},\omega_{0},q_{0},\tilde{\rho},(R_{a},D_{a})\}$ are the model parameters to be constrained, and ‘$pc$’ stands for the distance unit ‘parsec’. Here $(R_{a},D_{a})$ is the direction of anisotropy i.e., axis of cosmic rotation in equatorial coordinates. Further $m_{B,\rm corr}=m_{B}+\alpha x_{1}-\beta\texttt{c}+\gamma G_{\rm host}$ is the corrected apparent $B$-band magnitude of an SNIa from observations for its correlation with light curve fitting parameters viz., stretch ($x_{1}$), color (c) and host mass bias ($G_{\rm host}$). The corresponding coefficients $\alpha$, $\beta$ and $\gamma$ are global fit parameters along with the absolute $B$-band magnitude ‘$M_{0}$’ that are conventionally assumed to be same for all SNIa. ‘$M_{0}$’ is notably degenerate with the Hubble parameter $H_{0}$ (i.e., $h_{0}$), underscoring the importance of using SH0ES Cepheid variables in breaking that degeneracy.

The $\chi^{2}({\boldsymbol{\Xi}})$ function, which depends on various fitting parameters ${\boldsymbol{\Xi}}$, is defined as,

where the goal is to minimize this function to derive the desired constraints. Here, $\Delta{\boldsymbol{\mu}}=(\Delta\mu_{1},\Delta\mu_{2},\ldots,\Delta\mu_{n})^{T}$ is a column vector with ‘$n$’ elements corresponding to total number of SNe Ia in the dataset. The residual for the $i$-th SNIa is given by,

where,

and $\mu^{\rm th}_{i}$ denotes the theoretical distance modulus dependent on the model parameters for that $i$-th SNIa, . The above relation (Eq. 12) implies that whenever an SNIa host galaxy also contains a Cepheid, then the theoretical distance modulus $\mu^{\rm th}_{i}$ is replaced by the well measured value of distance modulus ‘$\mu^{\rm Cepheid}_{i}$’ to that galaxy obtained using the tight Cepheid period-luminosity relation. This way the degeneracy between the (dimensionless) Hubble parameter ‘$h_{0}$’ and the absolute magnitude ‘$M_{0}$’ is removed.

Apart from the global fit parameters ‘$\alpha$’ and ‘$\beta$’ corresponding to stretch and color correction of a supernova light curve, another dimensionless parameter ‘$\gamma$’ captures the dependence of apparent magnitude ‘$m_{B}$’ on host-galaxy stellar mass viz., $G_{{\rm host},i}$. This host-galaxy stellar mass parameter is defined as a step function depending on its mass as $G_{{\rm host},i}=+\frac{1}{2}$ if $M_{{\rm host},i}>10^{10}M_{\odot}$ and $G_{{\rm host},i}=-\frac{1}{2}$ if $M_{{\rm host},i}<10^{10}M_{\odot}$, where $M_{\odot}$ denotes 1 solar mass. The values of $x_{1,i}$, $\texttt{c}_{i}$, $m_{B}^{i}$, $G_{{\rm host},i}$ (as $\log_{10}(M_{{\rm host},i}/M_{\odot})$) and $\mu_{i}^{\rm Cepheid}$ are provided as part of Pantheon+ dataset for each SNIa.

Note that Pantheon+ collaboration also provides $m_{B,\rm corr}^{i}$ for each SNIa. Instead, we use $m_{B}^{i}$ and derive constraints on $\alpha$, $\beta$ and $\gamma$ along with cosmological parameters of our model. To account for the correlations among SNIa, which are influenced by both statistical and systematic uncertainties in the light curve fitting, we utilize the full covariance matrix. ${\bf C}_{\rm stat+sys}^{\rm SN}$ represents the covariance matrix for supernovae, while ${\bf C}_{\rm stat+sys}^{\rm Cepheid}$ denotes the covariance matrix for the SH0ES Cepheid-host distances, both of which are also provided as part of the Pantheon+ data release.

This section presents our results from an MCMC likelihood analysis of the Gödel-type rotating cosmological model, discussed in Sec. 2, to probe cosmic rotation. The full set of parameters to be constrained are ${\boldsymbol{\Xi}}=\{{\boldsymbol{\xi}},M_{0},\alpha,\beta,\gamma\}=\{h_{0},\Omega_{0},q_{0},\tilde{\rho},(R_{a},D_{a}),M_{0},\alpha,\beta,\gamma\}$ along with the cosmological parameters ${\boldsymbol{\Xi}}$ that describe the model of our present study. Since the luminosity-distance vs redshift relation given in Eq. (9) is a Taylor expanded relation up to linear order in redshift, $Z$, we derive cosmographic constraints on these parameters using low redshift SNIa data up to a maximum of $Z=0.5$. Specifically we estimate these parameters in five redshift bins ranging from $Z\leq 0.1$, $Z\leq 0.2$, $Z\leq 0.3$, $Z\leq 0.4$, and $Z\leq 0.5$. The astrophysical parameters $\{M_{0},\alpha,\beta,\gamma\}$ correspond to the SNIa absolute magnitude and coefficients of SNIa light-curve fitting parameters related to stretch, color, and host galaxy properties respectively. Actual cosmological parameters of our model are ${\boldsymbol{\xi}}=\{h_{0},\Omega_{0},q_{0},\tilde{\rho},(R_{a},D_{a})\}$. All these parameters are constrained using the latest Pantheon+ Type Ia supernova dataset.

In Fig. [2], we present the 2D contour plots illustrating the $1\sigma$ and $2\sigma$ confidence levels (CL) for the model parameters $h_{0},\Omega_{0},q_{0}$, and $\tilde{\rho}$ from different redshift bins mentioned above. It is important to note that these values should not be directly compared with $\Lambda$CDM constraints from SH0ES and Pantheon+ collaborations [69, 70], as the latter were obtained using the full dataset containing SNe Ia upto $Z\lesssim 2.3$. Thus, we also perform a similar cosmographic analysis of standard $\Lambda$CDM model. These results are presented in Fig. [3] to make appropriate comparison by selected red-shift bins. All the parameter constraints derived from the two models for the chosen redshift bins are tabulated in Table 1.

The Hubble constant ‘$h_{0}$’ in the Gödel-type model as well as $\Lambda$CDM model remains essentially same across all five redshift bins, with $h_{0}\approx 0.73$ for $Z\leq 0.1$, $\leq 0.2$, $\ldots$, $\leq 0.5$. These results agree well with the SH0ES determination of $H_{0}=73.04\pm 1.04$ (km/s/Mpc) [69]. The agreement suggests two possibilities, that the determination of ‘$h_{0}$’ using Pantheon+ SNe Ia compilation is driven by the tighter constraining power of SH0ES subset (that helps breaking the degeneracy between $h_{0}$ and $M_{0}$) or that the introduction of rotation in the Gödel-type framework may not be affecting the (local) expansion rate significantly.

The maximum likelihood values of the deceleration parameter ‘$q_{0}$’, shown in Fig. [2], for the Gödel-type model are $q_{0}\approx 0.01$, $-0.08$, $-0.077$, $-0.089$, and $-0.064$ for the five redshift ranges $Z\leq 0.1$, $\ldots$, $\leq 0.5$, respectively. A small positive value for $q_{0}$ at $Z\leq 0.1$ may be due to the limited constraining power of SNe Ia at low redshift or the influence of mild rotational effects on $q_{0}$ or nature of local matter distribution. It is clear though, from Table 1, that $q_{0}$ is consistent with ‘zero’ deceleration/acceleration when using SNe Ia upto $Z\leq 0.2$, given the error bars. However, as we include more number of supernovae with higher redshifts, it tends to take negative values consistent with $\Lambda$CDM predictions. Here we note that SH0ES team reported $q_{0}=-0.510\pm 0.024$ that includes some intermediate redshift SNe also (up to $Z<0.8$) in their data from the Pantheon+ sample [69]. In comparison, as shown in Fig. [3], the $\Lambda$CDM model with no rotation exhibits consistently negative values, inline with the standard model expectations. We find $q_{0}\approx-0.07$, $-0.22$, $-0.151$, $-0.117$, and $-0.084$ for the same redshift ranges, though $q_{0}$ is consistent with ‘zero’ given the error bars in the first redshift bin. We observe that these estimates of $q_{0}$ from $\Lambda$CDM model are smaller (more negative) than the values found with Gödel-type model for our Universe. Thus, this may be indicative of the (local) effects of rotation vis-a-vis expansion.

A key parameter in this study corresponding to global rotation, $\Omega_{0}$, exhibits a decreasing trend with increasing redshift bin. We obtain $\Omega_{0}=0.24\pm 0.46$ for $Z\leq 0.1$ that essentially implies no rotation, that then peaks to $\Omega_{0}=0.29^{+0.21}_{-0.15}$ for $Z\leq 0.2$, and thereafter reducing to $\Omega_{0}=0.182\pm 0.095$, $0.168\pm 0.072$, and $0.120^{+0.078}_{-0.046}$ for $Z\leq 0.3$, $Z\leq 0.4$, and $Z\leq 0.5$, respectively. The trend suggests that the effects of rotation, while uncertain at low redshifts, become more stringent as higher redshift data are included. Another dimensionless parameter ‘$\tilde{\rho}$’, defined as $\tilde{\rho}=\sqrt{\sigma/(k+\sigma)}$, that is related to global rotation (see Eq. 3) shows a similar decreasing trend. Constraints derived on both $\Omega_{0}$ and $\tilde{\rho}$ are presented in Table 1 for different redshift ranges considered. The higher the redshift bin, the smaller the value of $\tilde{\rho}$ assumes. Its values at posterior maxima varies from $\tilde{\rho}=0.0181$ to $0.0098$ from low redshift bin $Z\leq 0.1$ to the highest bin $Z\leq 0.5$.

This may be indicating two possibilities. On one hand it may be suggesting that we have to go beyond the cosmographic analysis. On the other hand, rotation due to local anisotropic nature of structure could be weakening as we go farther into the Hubble flow. A fuller study of this Gödel-type model, solving the formal luminosity vs redshift relation for any redshift, could provide a resolution to this ambiguity.

The constraints on absolute SNIa magnitude ‘$M_{0}$’ remain consistent being $\approx-19.2$ from the Gödel-type model as well as the $\Lambda$CDM model across all redshift ranges. Its value ranges from $M_{0}\approx-19.203$ to $-19.194$. This also validates the standard candle nature of SNe Ia. Here also it is important to note that we should not compare these values directly with the Panthoen+SH0ES determination of $M_{0}=-19.253\pm 0.027$ [69] that is estimated from the full data set. The astrophysical nuisance parameters $\alpha$, $\beta$, and $\gamma$ remain almost same across all redshift bins, reinforcing the robustness of SNIa light-curve corrections. Posterior distribution of these nuisance parameters along with $\tilde{\rho}$ and $M_{0}$ from our MCMC analysis are shown in Fig. [4].

Finally, we plot a derived parameter $\Omega_{0}/h_{0}$ ($=\omega_{0}/H_{0}$), as a 1D contour in Fig. 4 (bottom right panel), that captures the strength of rotation over expansion. The influence of rotation, though modest at $Z\leq 0.1$, becomes more pronounced for $Z\leq 0.2$ case where $\Omega_{0}/h_{0}\approx 0.29$. Thereafter rotation to expansion ratio decrease as $\Omega_{0}/h_{0}\approx 0.182,0.168$ and $0.120$ for $Z<0.3,0.4$, and $0.5$ respectively, indicating a diminishing impact of rotation as more distant SNe are included. This suggests that cosmic rotation, if present, is more relevant at lower redshifts but becomes negligible on larger scales. Here we note that, since $h_{0}$ is found to be nearly same ($\sim 0.73$) across all redshift ranges, the values of $\Omega_{0}/h_{0}$ show the same trend as $\Omega_{0}$ parameter.

In the end, we present one of the important results of this work viz., constraints on the global axis of anisotropy in our assumed rotating cosmological model. The inferred rotation axes are : $(R_{a},D_{a})\approx(251^{\circ},-41^{\circ})$ for $Z\leq 0.1$, $(244.8^{\circ},-46^{\circ})$ for $Z\leq 0.2$, $(239.6^{\circ},-56^{\circ})$ for $Z\leq 0.3$, $(239.3^{\circ},-52^{\circ})$ for $Z\leq 0.4$, and $(238.1^{\circ},-52^{\circ})$ for $Z\leq 0.5$ in equatorial coordinates that are also listed in Table 1. To put these anisotropy axes constraints in a broader perspective, these axes from different redshift bins with $1\sigma$ confidence contours are shown in galactic coordinates along with other anisotropy axes seen in various astronomical and cosmological data sets.

One of the earliest observational constraints on cosmic rotation was reported by Birch in 1982 where a dipole anisotropy in the distribution of observed angle between polarization and position angles of 137 radio sources was interpreted as evidence for cosmic rotation [31]. Birch reported the cosmic rotation axis as $(R_{a},D_{a})\approx(200^{\circ},-40^{\circ})$ in equatorial coordinates. Later, the same was confirmed using ‘indirectional statistics’ by Kendall and Young [76] with a higher significance who reported the axis of rotation to be $(R_{a},D_{a})\approx(202^{\circ},-37^{\circ})$ (using a sample of 134 sources) which is essentially same as that reported by Birch. Using a different sample of 205 sources, Andreasyan reported a cosmic anisotropy in the direction $(R_{a},D_{a})\approx(188^{\circ},-11^{\circ})$ [77]. Considering a subsample of Birch’s data (with 51 sources whose redshifts were known), a global axis of rotation within the present framework of the Gödel-like model was found to be $(R_{a},D_{a})\approx(184^{\circ},-38^{\circ})$ [30].

The present model of our interest given by Eq. (2) was earlier analyzed in Ref. [78] using the gold and silver Type Ia supernovae data sets of Ref. [79]. However, in that study, the rotation term (with the factor $\Omega_{0}/h_{0}$) in Eq. (9) was ignored constraining only the anisotropy axis (of rotation) and the parameter $\tilde{\rho}$. The gold data set had 157 SNIa and silver data set contained 29 sources. Because of the small sample size, perhaps, they could only find an upper bound on $\tilde{\rho}$ as $<0.049$ and $<0.046$, and the axis of anisotropy was found to be $(R_{a},D_{a})\approx(229.5^{\circ},-53.2^{\circ})$ and $(199.6^{\circ},-47.5^{\circ})$ in equatorial coordinates, when using the gold and gold+silver data sets respectively.

Future observations, particularly with additional intermediate and low redshift supernovae, will help refine these constraints and further assess the robustness of the anisotropy signal.

Rotation is a generic feature arising in any astrophysical setting that is an observational fact as seen in a wide variety of celestial systems. In this study, we explored the possibility of cosmic rotation for our Universe using a Gödel-type metric and observational data from Type Ia supernovae across multiple redshift ranges. We performed a cosmographic analysis of the magnitude vs redshift ($m-Z$) relation derived using the Kristian-Sachs formalism. We constrained both cosmological and astrophysical parameters using an MCMC method. Our results indicate interesting deviations from standard $\Lambda$CDM constraints for various redshift bins, particularly in the derived parameter $\Omega_{0}/h_{0}(=\omega_{0}/H_{0})$, that captures the strength of rotation vs expansion.

The inferred cosmological parameters, particularly the (dimensionless) Hubble constant $h_{0}$ and the deceleration parameter $q_{0}$, are central to understanding the effects of rotation in our analysis. We find that $h_{0}\sim 0.73$ being virtually same for all choices of redshift bins, indicating that the inclusion of rotational degrees of freedom in the Gödel-type model does not significantly alter the overall expansion rate compared to $\Lambda$CDM. However, the deceleration parameter exhibits notable variations across redshift bins, with values ranging from $q_{0}\approx 0.01$ at $Z\leq 0.1$ to $q_{0}\approx-0.089$ at $Z\leq 0.4$. In contrast, $q_{0}=-0.51$ for the standard flat $\Lambda$CDM model from Pantheon+ SNe Ia constraints [69]. The small positive value at very low redshift and the increasingly negative values at higher redshifts may arise due to statistical limitations in the dataset or subtle rotational effects influencing local expansion. As more high-redshift supernovae are included, $q_{0}$ approaches the standard $\Lambda$CDM expectation, albeit with slightly less negative values than those obtained in a purely isotropic model. This suggests that while rotation does not drastically impact cosmic acceleration, it may introduce small but noticeable deviations from the standard paradigm.

To quantify the relative impact of rotation versus expansion, we analyzed the dimensionless parameter $\Omega_{0}/h_{0}$, which provides insight into the strength of cosmic rotation across different redshift bins. Our results indicate that rotational effects are most pronounced at $Z\geq 0.1$, reaching a maximum value of $\approx 0.40$ for the $Z\leq 0.2$ bin. At lower redshifts ($Z<0.1$), this ratio is relatively small (essentially consistent with zero within error bars). However, as we include higher redshifts, the influence of rotation is found to diminish steadily, with $\Omega_{0}/h_{0}$ decreasing to $\approx 0.16$ for $Z\leq 0.5$ bin. This decline suggests that rotational effects, if present, are confined to relatively local cosmic scales and become negligible when considering more distant supernovae. These results highlight that while rotation may play a role in the cosmic evolution at small redshifts, its contribution becomes increasingly negligible on larger scales, reinforcing the overall isotropy of the universe. The global rotation axis inferred from our model is consistent across all redshift bins, with its mean direction given by $(R_{a},D_{a})\approx(243^{\circ},-49^{\circ})$.

A related parameter that captures global rotation in this model is $\tilde{\rho}=\sqrt{\sigma/(k+\sigma)}$. It follows the same trend as $q_{0}$ and $\Omega_{0}$ (and hence $\Omega_{0}/h_{0}$). It values decreases steadily from the lowest ($Z\leq 0.1$) to highest ($Z\leq 0.5$) redshift ranges chosen in the present study. Astrophysical parameters viz., $M_{0}$, $\alpha$, $\beta$, and $\gamma$, which characterize supernova magnitude calibration, its light-curve stretch, color and host galaxy properties, are nearly same across different redshift ranges. This reinforces the robustness of the dataset and suggests that systematic effects from supernova standardization are minimal.

The AIC-based model comparison reveals an intriguing preference for the Gödel-type rotating model over the standard flat $\Lambda$CDM cosmology, particularly at higher redshifts. While both models remain competitive for $Z\leq 0.1$, with the candidate model being only marginally more probable, the preference for the rotating model increases significantly at larger redshift cutoffs. Specifically, for $Z\leq 0.3$, $0.4$, and $0.5$, the Gödel-type model is found to be $\sim 9.9$, $4.7$, and $6.9$ times more probable than $\Lambda$CDM model, respectively. This suggests that the additional kinematic effects introduced by cosmic rotation might play a non-negligible role in the late-time expansion history. However, caution is warranted in interpreting these results. The robustness of the cosmographic approach and potential systematics unaccounted for in the data could have a bearing on these findings and further investigation - fully testing this model going beyond the cosmographic expansion - are necessary to determine whether this preference is indicative of a true physical effect or a statistical fluke or a data artifact.

In conclusion, our analysis shows indications for cosmic rotation at intermediate reshifts given that the AIC model assessment metric strongly favours the Gödel-type model over the standard model. Nevertheless it needs further study to derive better constraints on the model to fully understand the ramifications of cosmic rotation more clearly. Future observations from upcoming surveys like Euclid [82] and JWST [83] may provide improved constraints on both the rotation parameter and the anisotropy axis. This could offer new insights into cosmic isotropy and potential deviations that enrich our understanding of the universe’s large-scale structure.

AV acknowledges the financial support received through a research fellowship awarded by Council of Scientific & Industrial Research (CSIR), India during this project. We acknowledge National Supercomputing Mission (NSM) for providing computing resources of ‘PARAM Shivay’ at Indian Institute of Technology (BHU), Varanasi, which is implemented by C-DAC and supported by the Ministry of Electronics and Information Technology (MeitY) and Department of Science and Technology (DST), Government of India. DFM thanks the Research Council of Norway for their support and the resources provided by UNINETT Sigma2 – the National Infrastructure for High-Performance Computing and Data Storage in Norway. We acknowledge using Cobaya²²2https://cobaya.readthedocs.io/en/latest/ [84] - a code for Bayesian analysis in Cosmology, and GetDist³³3https://getdist.readthedocs.io/en/latest/ [85] - a user friendly GUI for the analysis and plotting of MCMC samples. We also acknowledge the use of SciPy⁴⁴4https://scipy.org/ [86], NumPy⁵⁵5https://numpy.org/ [87], Astropy⁶⁶6http://www.astropy.org [88, 89, 90], and Matplotlib⁷⁷7https://matplotlib.org/ [91].