An Accurate Fitting Function For Scale-dependent Growth Rate in Hu-Sawicki $f(R)$ Gravity

TL;DR

This paper introduces a highly accurate, simple fitting formula for the scale-dependent growth rate in Hu-Sawicki $f(R)$ gravity, enabling precise predictions for cosmological observations across various scales and redshifts.

Contribution

The paper provides a new fitting function for the growth rate in Hu-Sawicki $f(R)$ gravity with sub-percent accuracy, validated over a wide range of scales, redshifts, and model parameters.

Findings

Achieves sub-percent maximum error in fitting the growth rate.

Valid across scales $10^{-4} \,\leq k \leq 5$ Mpc$^{-1}$ and redshifts $0 \leq z \leq 3$.

Useful for cosmological probes like redshift space distortions and galaxy velocities.

Abstract

We present a simple fitting formula for the scale-dependent growth rate of Hu-Sawicki model in $f(R)$ modified gravity. We compare the accuracy of the fitting function against numerical results and report achieving a sub-percent maximum error for all different studied cases. The validity of the fitting function is tested against a wide range of scales, $10^{-4} \le k \le 5 \ [\mathrm{Mpc}^{-1}]$ , for various redshift values in the range of $z \in [0.0,\ 3.0]$ and different $|f_R^0|$ values. This formula is useful in producing predictions for cosmological probes that are sensitive to growth rate such as redshift space distortions, galaxy-galaxy correlation function, peculiar velocity statistics and pairwise velocities of galaxy clusters.

Equations

$$S=\int d^{4}x\sqrt{-g}\bigg{[}\frac{f(R)}{16\pi G}+\mathcal{L}^{(m)}\bigg{]},$$

$$F{R}_{\mu \nu }-\frac{1}{2}f{g}_{\mu \nu }+(g_{\mu \nu }\square -{\nabla }_{\mu }{\nabla }_{\nu })F=8\pi G{T}_{\mu \nu },$$

$$F(R)\equiv \frac{\partial f(R)}{\partial R}.$$

$$M^{2}=\frac{1}{3}\bigg{(}\frac{F}{F_{,R}}\bigg{)},$$

$${\lambda }_c=a^{-1}\sqrt{3\frac{\partial F(R)}{\partial R}}.$$

$$\ddot{\delta}_{m}+\bigg{(}2H+\frac{\dot{F}}{2F}\bigg{)}\dot{\delta}_{m}-\frac{\rho_{m}}{2F}\delta_{m}=\frac{1}{2F}\bigg{[}\Big{(}-6H^{2}+\frac{k^{2}}{a^{2}}\Big{)}\delta F+3H\dot{\delta F}+3\ddot{\delta F}\bigg{]},$$

$$\ddot{\delta F}+3H\dot{\delta F}+\bigg{(}\frac{k^{2}}{a^{2}}+\frac{F}{3F_{,R}}-\frac{R}{3}\bigg{)}\delta F=\frac{1}{3}\rho_{m}\delta_{m}+\dot{F}\dot{\delta}_{m}\ .$$

$${\ddot{\delta }}_m+2H{\dot{\delta }}_m-4\pi G_{\mathrm{eff}}{\rho }_m{\delta }_m\simeq 0,$$

$$G_{\mathrm{eff}}(k,z)=\frac{G}{F}\bigg{[}1+\frac{1}{3}\Big{(}\frac{k^{2}}{a^{2}M^{2}/F+k^{2}}\Big{)}\bigg{]}.$$

$$f(R)=R-2\Lambda \frac{R}{R+{\mu }^2},$$

$$f(R)=R-2\Lambda -f_R^0\frac{{\bar{R}}_0^2}{R},$$

$$\lambda_{C}=(1+z)\sqrt{6|f_{R}^{0}|\bigg{(}\frac{\bar{R}_{0}^{2}}{R^{3}}\bigg{)}}.$$

$$f_g(a)=\frac{d\ln {\delta }_m(a)}{d\ln a}.$$

$$f_g(a)={\Omega }_m^{\gamma }(a).$$

$$f_g(k,a)=\frac{d\ln {\delta }_m(k,a)}{d\ln a}.$$

$$f_{g}^{\prime}(k,z)-\frac{f_{g}^{2}(k,z)}{1+z}+\Big{(}\frac{E^{\prime}(z)}{E(z)}-\frac{2}{1+z}\Big{)}f_{g}(k,z)+\frac{3}{2}\Omega_{m}(z)(1+z)\frac{G_{\mathrm{eff}}(k,z)}{G}=0,$$

$$\large f_{g}(k,z;\ |f_{R}^{0}|)=\Omega_{m}^{\gamma}(z)\bigg{(}1+\frac{e^{\frac{-\alpha(z)}{\mathcal{K}}}}{5.5}\bigg{)},$$

$$\mathcal{K}\equiv k{\lambda }_C.$$

$$\alpha (z)$$

$$\Delta f_g(\%)\equiv (\frac{f_g^{fit.func.}-f_g^{\mathrm{numerical}}}{f_g^{\mathrm{numerical}}})\times 100.$$

$$P_s(k,\mu )=(1+\beta {\mu }^2{)}^2{P}_r(k).$$

$$\beta (k,z)=\frac{f_g(k,z)}{b(z)},$$

$$v(\mathbb{k})=-i{H}_0{f}_g\frac{\hat{\mathbb{k}}}{k}{\delta }_m(\mathbb{k}).$$

$${\Psi }_{ij}(r)\equiv ⟨v_i(x)v_j(x+r)⟩,$$

$$\Psi_{ij}(r)=\Psi_{\bot}(r)\delta_{ij}+\big{[}\Psi_{||}(r)-\Psi_{\bot}(r)\big{]}\hat{r}_{i}\hat{r}_{j}$$

$${\Psi }_{\perp }(r)=\frac{H_0^2}{2{\pi }^2}\int dk{f}_g^2(k,z)P(k)\frac{j_1(kr)}{kr},$$

$${\Psi}_{||}(r)=\frac{H_{0}^{2}}{2\pi^{2}}\int dkf_{g}^{2}(k,z)P(k)\bigg{[}j_{0}(kr)-2\frac{j_{1}(kr)}{kr}\bigg{]},$$

$$C^{g_i{g}_j}(\ell )=\int dz\frac{H(z)}{{\chi }^2(z)}{b}_i{W}_{g_i}(z)b_j{W}_{g_j}(z)[D_m^2(k,z)P(k){]}_{k=\frac{\ell +1/2}{\chi (z)}}.$$

$$\int d\ln D_m(k,z)=-\int dz\frac{f_g(k,z)}{1+z},$$

Full text

1 Introduction

Since its unexpected discovery more than 20 years ago [36, 33] the accelerated expansion of the Universe has remained in the spotlight as one of the major unsolved problems of modern physics. Cosmological measurements of distant type Ia supernovae were the original indicators of the existence of this phenomenon. Subsequently, these discoveries were followed by observations in cosmic microwave background anisotropies and large scale structure surveys, all supporting the aforementioned discovery [2]. The standard paradigm of modern cosmology, $\Lambda\mathrm{CDM}$ model, explains the late time accelerating expansion of the Universe by virtue of an additional constant, $\Lambda$, added to the Einstein field equations. However, the discrepancy between the observed and the theoretically expected value for the cosmological constant leads to different phenomenological and philosophical problems [6, 41]. Hence, several theoretical approaches and new frameworks have been proposed to solve this problem [19]. The majority of the developed solutions belongs to one of the two different paradigms: 1 - Dark Energy models, where exotic new forms of energy are introduced to source the late time accelerated expansion of the Universe. The cosmological constant can be classified as a specific member of this group. 2 - Modified Gravity models, that are based on the assumption of the inadequacy of General Relativity (GR) to explain the dynamics of the gravitational field on certain cosmological scales. Here, the theory of gravity changes and typically a new degree of of freedom is introduced in the theory that plays the role of the cosmological constant naturally [9].

The distant supernovae observations and the background expansion history alone are not sufficient to distinguish between different Dark Energy and Modified Gravity models [37]. On the other hand the growth of cosmological perturbations has been shown to have the potential to discriminate between the two models [34, 26]. Upcoming surveys like Euclid, DESI, WFIRST, LSST [3, 1, 16] will map huge volumes of the Universe and the structure within, opening the door for detailed analysis of the structure formation over the last three quarters of the history of the Universe [15] which in turn will give the researchers the opportunity to better discriminate between alternative models to cosmological constant scenario.

Motivated by these arguments, we present a fitting formula for the scale-dependent growth rate of $f(R)$ gravity. The fitting function is capable of reproducing the numerical results in redshifts and scales that will be measured accurately by the upcoming and future surveys. Previous attempts to model the growth rate of modified gravity with analytical or semi-analytical approaches have been able to reproduce acceptable results [40]. Here we present a new fitting function that is capable of reproducing the numerical results up to very small scales, namely $10^{-4}\leq k\leq 5\ [\mathrm{Mpc}^{-1}]$ for the redshift range of $z\in[0.0,3.0]$ with a sub-percent accuracy. In addition, the explicit dependence of the fitting formula on $|f_{R}^{0}|$ provides a simple but accurate way of using the growth rate in future cosmological analysis.

This paper is organized as follows. In section 2 we briefly review the $f(R)$ gravity and the specific model that we have used in this study. Section 3 is a short review of structure formation for $\Lambda\mathrm{CDM}$ and $f(R)$ gravity with the emphasise being on the growth rate. In section 4 we present the fitting formula for the growth rate of matter density perturbations and scrutinize it against different cosmological scenarios. Section 5 is dedicated to potential applications of the fitting function in large-scale structure cosmology. Lastly, Section 6 is our final remarks and conclusions. Throughout the paper we have used $H_{0}=70\ [\mathrm{km\ s^{-1}Mpc^{-1}}]$ while assuming a flat Friedmann-Lemaitre-Robertson-Walker (FLRW) background.

2 $f(R)$ Gravity Review

A potential model for explaining the late time accelerated expansion of the Universe is modified gravity. Among different proposed models, $f(R)$ modified gravity is of particular interest. By employing the chameleon mechanism [27, 22] these theories modify General Relativity on cosmological scales while at the same time passing the local tests of gravity. In the $f(R)$ model of gravity the Einstein-Hilbert action is replaced by

[TABLE]

where $G$ is the Newton’s gravitational constant and $\mathcal{L}^{(m)}$ is the matter Lagrangian [30, 39]. $f(R)$ is a function of the Ricci scalar, $R$ and the specific functional form is determined by a combination of theoretical and phenomenological considerations. In the special case of $f(R)=R-2\Lambda$ one recovers the Einstein-Hilbert action plus cosmological constant and consequently $\Lambda\mathrm{CDM}$ cosmology. Modified Einstein field equations have the following form in $f(R)$ gravity:

[TABLE]

where

[TABLE]

$F$ is a new degree of freedom of the theory that acts similar to a scalar field, dubbed scalaron, with a mass given by

[TABLE]

where $F_{,R}=\partial^{2}f(R)/\partial R^{2}$. The scalaron field propagates on scales that are smaller than its associated comoving Compton wavelength,

[TABLE]

Modified gravity models in general are able to reproduce and mimic the redshift-distance measurements of supernove and are indistinguishable from $\Lambda\mathrm{CDM}$ and Dark Energy models in this regard [4]. Therefore, "background"-level observations alone are not sufficient to test or falsify these models. Despite that, in linear perturbation regime the equations governing the relationship between matter overdensities and the gravitational potential are also modified, leading to a different growth characteristic and growth history. Hence modified gravity have unique and falsifiable predictions in this context and can be tested against future observations. In particular, $f(R)$ models modify the growth of density perturbations and structure formation via introduction of a scale-dependant growth which is in contrast to $\Lambda\mathrm{CDM}$ and the scale-independent growth therein [40].

At linear perturbation level, in the comoving gauge the modified Einstein equations lead to the following perturbation equations in Fourier space for the evolution of matter overdensities:

[TABLE]

[TABLE]

Here $k$ is the comoving wavenumber, $a=(1+z)^{-1}$ is the scale factor normalized to unity for the current time, $\rho_{m}$ is the matter density, $\delta_{m}(a)=\delta\rho_{m}/\rho_{m}$ is the density contrast, $H=\dot{a}/a$ is the Hubble parameter and "\ {}^{.}\" represents derivative with respect to the cosmic time. Since we are interested in the evolution of density perturbations in $f(R)$ gravity we combine eq. $eq.~{}\eqref{perturb1}$ and $eq.~{}\eqref{perturb2}$ to obtain the modified equation of evolution for the matter density contrast. For $f(R)$ models that are cosmologically viable, F has small variation ($|\dot{F}|\ll HF$), consequently in both equations $\dot{F}$ can be neglected. The oscillating mode is also negligible compared to modes induced by matter perturbations. Finally since we are interested in modes well within the Hubble radius we have $k^{2}/a^{2}\gg H^{2}$. All the above assumptions leads to the following equation,

[TABLE]

where $G_{\mathrm{eff}}$ is the "effective" gravitational constant, defined as

[TABLE]

Since $RF_{,R}\ll 1$, $eq.~{}\eqref{mass}$ can be rewritten as $M^{2}=F/3F_{,R}$. It is worth noting that $G_{\mathrm{eff}}$ has replaced $G$, which would instead appear in standard $\Lambda\mathrm{CDM}$ cosmology. Since $G_{\mathrm{eff}}$ is a scale-dependent function, the structure formation in modified gravity includes scale dependent features that are unique and can in principle be used to distinguish $f(R)$ models from standard cosmology [29].

So far we have not assumed any specific functional form for $f(R)$. A viable model should be able to satisfy a wide range of cosmological and local constraints imposed by observational and theoretical considerations. In this work we use the Hu-Sawicki $f(R)$ model (with $n=1$) that has been previously shown to pass the local solar system tests [14]. The functional form for this model is as follows,

[TABLE]

where $\Lambda$ and $\mu$ are free parameters. With the assumption of small deviations from GR ($R\gg\mu^{2}$), $eq.~{}\eqref{hu}$ can be simplified as:

[TABLE]

where $|f_{R}^{0}|=2\Lambda\mu^{2}/\bar{R}_{0}^{2}$. Here, $R=6(2H^{2}+\dot{H})$ and $\bar{R}_{0}=\bar{R}(z=0)$ where overbar indicates a background quantity. Using eq. $\eqref{eq:2.11}$, the comoving Compton wavelength for this model is

[TABLE]

It is conventional to use $|f_{R}^{0}|$ to quantify the deviation of $f(R)$ models from General Relativity. A smaller value of $|f_{R}^{0}|$ is an indication of weaker deviation. Latest cosmological studies indicate that $|f_{R}^{0}|\leq 10^{-5}$ [7, 25]. Regardless, we have assumed a wide range of potential values for $|f_{R}^{0}|=\{10^{-6},10^{-5},10^{-4}\}$ in this work. The inclusion of $|f_{R}^{0}|=10^{-4}$ in our analysis is for demonstration purposes only.

3 Growth Rate of Matter Density Perturbations

The growth rate of matter density perturbations in $\Lambda\mathrm{CDM}$ cosmology is a function of scale factor alone [32, 23] and is defined as,

[TABLE]

It has been shown that for $\Lambda\mathrm{CDM}$ cosmology, to a high level of accuracy, $f_{g}(a)$ can be approximated by

[TABLE]

where $\gamma\simeq 0.55$ is the growth index and $\Omega_{m}(a)=\Omega_{m,0}a^{-3}E^{-2}(a)$ is the matter density for flat FLRW cosmology [42] with $E(z)=\frac{H(z)}{H_{0}}=\sqrt{\Omega_{m,0}(1+z)^{3}+\Omega_{\Lambda}}$. In the case of $f(R)$ gravity, due to the scale-dependent nature of $G_{\mathrm{eff}}(k,z)$, the growth rate is expected to have scale-dependent features as well and can be written as,

[TABLE]

Since our aim is to study the growth rate of matter perturbations in $f(R)$ gravity, using $eq.~{}\eqref{fgfR}$ we can rewrite eq. $eq.~{}\eqref{eq:2.6}$ in terms of $f_{g}$ and in a form that is numerically more tractable,

[TABLE]

where we have changed the differentiation variable from cosmic time to redshift, ${}^{\prime}\equiv\partial/\partial z$. By solving $eq.~{}\eqref{eq:3.4}$ numerically, in Fig. 1 we show the scale-dependant behavior of growth rate in $f(R)$ gravity. The figure shows that deviations from $\Lambda\mathrm{CDM}$ become larger on small scales. In addition, a larger value of $|f_{R}^{0}|$ produces larger deviations. In Fig. 2 we plot the growth rate against redshift while having the scale fixed. As evident, the deviation from GR is a combination of redshift and scale.

4 Fitting Function For Growth Rate

Analytical and/or semi-analytical approaches to model the scale-dependent growth rate for modified gravity has been proposed before [35, 5, 29, 8]. We propose a simple fitting function that is capable of recreating the numerical results with high accuracy on a wide range of scales for various redshifts while having an easy to use functional form. We have calculated $f_{g}$ over a range of values for $|f_{R}^{0}|$ across different redshifts. Our intention is to find a fitting function that properly reproduce the numerical results not only on larger scales but also on small ones. Consequently, we present the following fitting formula

[TABLE]

where we have defined the new variable $\mathcal{K}$ as

[TABLE]

As before, $k$ is the comoving wave number and $\lambda_{C}$ is the comoving range of the scalar field force, defined in the previous section. $\alpha(z)$ is a redshift dependent functions defined as

[TABLE]

The truncation of $\alpha(z)$ at second order in terms of redshift is based on the observation that a linear truncation is not sufficient enough to produce the required accuracy (below 1$\%$ error), while at the same time, inclusion of a cubic term has no meaningful contribution to the results. We evaluated the accuracy of the fitting formula for different values of redshift, scale, $\Omega_{m}$ and $|f_{R}^{0}|$. Namely, we have considered cases where $z\in[0.0,3.0]$, $|f_{R}^{0}|\in\{10^{-6},10^{-5},10^{-4}\}$ and $10^{-4}\leq k\leq 5\ \mathrm{[Mpc^{-1}}]$ with $\Omega_{m}\in[0.28,0.335]$. We report that the fitting function is able to match the numerical results by a sub-percent accuracy for all the studied cases. Fig. 3 is a demonstration of the quality of the fitting function when compared against a large collection of numerical results. It is worth mentioning that the explicit dependence of the fitting formula on $|f_{R}^{0}|$ makes the computation of derivatives with respect to $|f_{R}^{0}|$ and consequently potential use in Fisher Matrix analysis straightforward. We assess accuracy of the fitting formula against numerical results by

[TABLE]

In Fig. 4, we plot the relative error of the fitting formula against numerical results. As evident in Fig. 4, the fitting formula has a sub-percent accuracy over a very wide range of scales for $z\in[0.0,\ 3.0]$ for different $|f_{R}^{0}|$ values. Here $\Omega_{m,0}=0.2815$ has been used.

5 Cosmological Probes and Applications

In this section we briefly discuss cosmological probes that are sensitive to the growth rate directly or indirectly where the fitting formula can be used to perform calculations for cosmological analysis.

5.1 Redshift Space Distortion

Clustering patterns of galaxies in redshift space is affected by their peculiar velocities. On large scales this phenomenon is called the Kaiser effect [21]: Infalling of galaxies into regions with high-density such as galaxy clusters. On small scales, the Fingers-of-God effect is present: The clustering pattern of galaxies is elongated along the line of sight in the redshift space, caused by random peculiar velocities of galaxies that are bound in virialised structures. In linear perturbation theory, the relation between the redshift space galaxy power spectrum and the power spectrum in real space is as follows [17],

[TABLE]

Here $\mu=\frac{\mathbb{k}\cdot\mathbb{r}}{kr}$ is the cosine of the angle between the line of sight and direction of infall. The distortion parameter, $\beta$, is defined as

[TABLE]

where $b(z)=\delta_{g}/\delta_{m}$ is the linear galaxy bias. In contrast with $\Lambda\mathrm{CDM}$, because of the scale dependent nature of growth in $f(R)$ gravity the distortion parameter is scale dependent too. In Fig. 5 we plot the distortion parameter against redshift for different $k$ modes for various $|f_{R}^{0}|$ values. We have assumed a simple and only redshift-dependent bias, $b(z)=\sqrt{1+z}.$

5.2 Velocity Correlation Function

In linear perturbation theory Fourier components of peculiar velocity and density contrast field are related as follows [31],

[TABLE]

In the context of $f(R)$ modified gravity since the energy-momentum conservation equations are unchanged, eq. (5.3) is not altered and can be used. In order to study the correlation of peculiar velocity field the velocity correlation tensor is defined as,

[TABLE]

where $\mathbb{r}$ is the separation vector between two points located at $r_{i}$ and $r_{j}$. For a statistically homogeneous and isotropic velocity field the correlation tensor is a linear combination of the parallel, $\Psi_{||}$ and transverse (to the separation vector), $\Psi_{\bot}$ components [20, 11],

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. It has been shown that given the matter power spectrum, $P(k)$, $\Psi_{||,\bot}$ can be calculated directly as

[TABLE]

and

[TABLE]

where $j_{0}(x)=\frac{sin(x)}{x}\hskip 1.42262pt$ and $j_{1}(x)=\frac{sin(x)}{x^{2}}-\frac{cos(x)}{x}$ are the spherical Bessel functions. We like to emphasize the fact that in standard $\Lambda\mathrm{CDM}$ case the scale-independent growth rate, $f_{g}(z)$ is not part of the integration. Since $\Psi_{||}$ and $\Psi_{\bot}$ are functions of power spectrum and growth rate, in principle one can study the parallel and transverse components of the velocity correlation tensor to constrain modified gravity. In Fig. 5 we have plotted $\Psi_{\bot}(r)$ and $\Psi_{||}(r)$ for $\Lambda\mathrm{CDM}$ and $f(R)$. As expected the $f(R)$ curves start deviating from GR on smaller scales. We have used the publicly available MGCAMB code [13, 43] to calculate the $f(R)$ power spectrum for different $|f_{R}^{0}|$ values.

5.3 Pairwise Velocity of Galaxy Clusters

The recent discovery of pairwise kinematic Sunyaev-Zel’dovich effect [12] is a new opportunity to use the second order CMB physics to study gravity. pkSZ can act as a proxy to study pairwise velocity of galaxy clusters and it has been shown to have great potential to constrain gravity and dark energy models on cosmic scales [28]. We are planning to address this problem in a more detailed fashion in an upcoming paper [10].

5.4 Galaxy-Galaxy Angular Correlation

By ignoring the very small angular multipoles, namely for $\ell\geq 10$ and using Limber approximation [24] the galaxy-galaxy angular auto-correlation can be simplified as,

[TABLE]

The growth function, $D_{m}(k,z)$, is related to the growth rate and can be expressed as

[TABLE]

where $f_{g}(k,z)$, the fitting formula can be used to perform the calculation. In equation (5.8), $\chi(z)$ is the comoving distance to redshift $z$, $b_{i}$ is the galaxy bias for the redshift bin $i$ and $W_{g_{i}}(z)$ is the galaxy redshift distribution in the same bin (normalized to 1). In addition to the galaxy-galaxy angular auto-correlation, the fitting formula can be used for calculation of angular cross-correlation between Cosmic Microwave Background temperature and galaxy distribution (Integrated Sachs-Wolfe cross-correlated with galaxies). Here, on top of using the fitting formula for obtaining $D_{m}(k,z)$, the integrated Sachs-Wolfe effect needs to be modeled in modified gravity context [38].

6 Conclusions

Future high precision large-scale structure surveys will map the distribution and evolution of matter densities in the Universe with unprecedented detail and will provide opportunities to study gravity and alternative models of gravity on cosmological scales. In this work we reviewed the scale-dependant growth of large-scale structure in $f(R)$ gravity under different assumptions about the underlying model and for various redshifts. Motivated by the potential of the future experiments to constrain modified gravity, we introduced a fitting formula for scale-dependent growth rate of matter perturbations that is in sub-percent agreement with numerical results over a wide range of scale, redshift and has explicit dependence on $|f_{R}^{0}|$, the parameter indicating deviation from $\Lambda\mathrm{CDM}$ cosmology. We discussed number of cosmological probes that directly or indirectly can be used to study the growth rate. Namely, the fitting formula is able to regenerate the numerical results with high accuracy for $10^{-4}\leq k\leq 5\ [\mathrm{Mpc}^{-1}]$, $0\leq z\leq 3$ for $|f_{R}^{0}|=\{10^{-6},10^{-5},10^{-4}\}$. Using the fitting formula, we calculated the peculiar velocity statistics for $f(R)$ gravity and compared the results to standard cosmology. Depending on the scale that is considered and the value for $|f_{R}^{0}|$ the $f(R)$ predictions can deviate from standard cosmology significantly. Lastly, we reviewed the redshift space distortion phenomena in the context of $f(R)$ gravity while using the fitting formula for obtaining the results.

Acknowledgments

We like to thank Shant Baghram and Siavash Yasini for helpful discussions. NM and EP are supported by NASA grant 80NSSC18K0403. NM was partially supported by the WiSE major support for faculty, awarded to EP. The authors acknowledge using the SciPy library [18] for preparation of the results presented in this work.