Axion-like particle (ALP) portal freeze-in dark matter confronting ALP search experiments

The relic density of Dark Matter (DM) in the freeze-in scenario is highly dependent on the evolution history of the universe and changes significantly in a non-standard (NS) cosmological framework prior to Big Bang Nucleosynthesis (BBN). In this scenario, an additional species dominates the energy budget of the universe at early times (before BBN), resulting in a larger cosmological expansion rate at a given temperature compared to the standard radiation-dominated (RD) universe. To investigate the production of DM in the freeze-in scenario, we consider both standard RD and NS cosmological picture before BBN and perform a comparative analysis. We extend the Standard Model (SM) particle content with a SM singlet DM particle $\chi $ and an axion-like particle (ALP) $a$. The interactions between ALP, SM particles, and DM are generated by higher dimensional effective operators. This setup allows the production of DM $\chi$ from SM bath through the mediation of ALP, via ALP-portal processes. These interactions involve non-renormalizable operators, leading to ultraviolet (UV) freeze-in, which depends on the reheating temperature ($T_{RH}$) of the early universe. In the NS cosmological scenario, the faster expansion rate suppresses the DM production processes, allowing for enhanced effective couplings between the visible and dark sectors to satisfy the observed DM abundance compared to RD scenario. This improved coupling increases the detection prospects for freeze-in DM via the ALP-portal, which is otherwise challenging to detect in RD universe due to small couplings involved. Using an effective field theory set-up, we show that various ALP searches such as in FASER, DUNE, and SHiP, etc. will be able to probe significant parameter space depending on the different model parameters.


I. INTRODUCTION
A variety of independent astrophysical observations have established the presence of an ample non-baryonic as well as electrically neutral dark matter (DM) constituent in our universe [1][2][3][4]. Additionally, from cosmological point of view, the presence of DM component of the universe is expected to play a significant role in the large scale structure formation of our universe [4]. However, the origin and composition of the particle nature of DM hypothesis is a modern day intriguing puzzle for both theoretical and experimental particle physics world so to say. Since the standard model (SM) with its current particle contents can not accommodate DM, several well motivated beyond the SM (BSM) scenarios suggest a suitable candidate for DM [5][6][7]. One such popular candidate is the weakly interacting massive particle (WIMP), which is attracts a lot of attentions in the context of thermal dark matter models [8,9]. The WIMP scenario demands dark sector-SM interaction cross-section of the order of electro-weak (EW) interaction strength (also known as WIMP Miracle) [10][11][12]. For the past several decades multiple DM search experiments have provided null results, thus pushing a very stringent bounds on the interaction strength between DM and the SM particles from direct detection experiments [13][14][15][16][17][18], indirect detection [19,20] and through collider searches, as in for example, at the Large Hadron Collider (LHC) [21][22][23].
In order to avoid these bounds on the WIMP scenario, alternative avenues where DM abundance may have been created out of equilibrium by the so-called freeze-in mechanism [24][25][26]. In this framework, DM particle is a feebly interacting massive particle (FIMP) having extremely suppressed coupling O(10 −12 − 10 −10 ) with the SM particles. Thus in this scenario, if DM exists it resides in a dark sector with such a tiny coupling with the SM particles so that DM particles never maintains thermal equilibrium with SM bath in the early Universe. Rather the DM particles get produced non-thermally from the decay or via annihilation of the SM bath particles. Thus from an initial negligible number density DM abundance increases and freezes in after sometime resulting a fixed relic abundance. This freeze in scenario is broadly classified into (a) Infra-red (IR) freeze-in relevant at lower temperature [24][25][26][27][28][29][30][31][32] and (b) Ultra-violate (UV) freeze-in which occurs at higher temperatures (approximately the reheating temperature of the Universe) [33][34][35][36][37][38][39][40][41]. The main challenging aspect of this scenario to see any laboratory imprints (direct searches or otherwise) due to the smallness of the coupling strengths (therefore weakly interacting dark sector) involved in the processes. However, some possibilities have been proposed to probe such scenarios. For example, if the DM production in the early Universe proceeds via the decay or scattering of thermal bath particles then feeble couplings associated with such decays would make the dark sector particles very long-lived and can be looked at the LHC and beyond [42][43][44][45][46][47][48][49][50][51][52]. Recently Gravitational Waves from early universe has also been proposed to test such freeze-in DM scenarios and weakly coupled dark sector in general [53,54].
Taking a slight detour from the ongoing discussion let us mention another shortcomings of the SM, which involves the strong CP problem which states that the θ-parameter of QCD is highly constrained (| θ | 10 −10 ) from measurements of the neutron electric dipole moment. It has been postulated by Peccei and Quinn (PQ) that to resolve this puzzle, one has to introduce a hypothetical particle known as axion [55][56][57][58] which is a (pseudo) Nambu-Goldstone boson of a spontaneously broken global chiral U (1) symmetry (known as PQ symmetry and denoted by U (1) PQ ) at some high scale f a >> Λ QCD ≈ O(100) MeV. On the other hand, the QCD anomaly explicitly breakes this U (1) PQ at the quantum level, leading to a tiny QCD axion mass, m a ∼ m π f π /f a , where m π is the pion mass and f π is the pion decay constant. In general, QCD axion mass and the magnitude of its couplings to the SM particles are inversely proportional to the axion decay constant f a . The main characteristic features of the QCD axion, like smallness of its mass and weak derivative couplings with the SM particles arise from its pseudo-Nambu Goldstone property. However, exact couplings of the QCD axion to the SM particles are highly model dependent and there are three well known UV complete models of QCD axion, namely PQWW [55][56][57], KSVZ [59,60] and DFSZ [61,62] have been extensively investigated in the literature. Among these, the PQWW axion model, where f a ∼ v SM , where v SM = 246 GeV is the SM vev, and m a ∼ 10 KeV is already ruled out from experimental measurement of rare karon decays, K + → π + a [63]. On the other hand for large f a >> v SM the axion is very light and very weakly coupled with the SM sector and can be long lived too. These are invisible axions and can be a suitable candidate for dark matter (DM) satisfying the observed relic density of DM via the misalignment mechanism [64][65][66][67][68][69]. For example, with f a ≈ 10 11 GeV and m a = 6 × 10 −5 eV, axions can accommodate ∼ 1% of the observed DM relic density [70]. Pseudo-Nambu Goldstone bosons (pNGBs) may appear in phethora of beyond the standard model scenarios, typically those having spontaneously broken one or more global symmetries at some high scale. Such pNGBs are broadly called axion like particles (ALP) . Unlike QCD axions, ALPs are neither required to solve the strong CP problem nor their mass and couplings are related to the QCD anomaly. Consequently, one can treat masses and couplings of ALPs as independent quantities, thus allowing for larger parameter space that can be probed and more scope for the study of model realizations. Like invisible QCD axion, ALPs are also valid dark matter candidates [71][72][73][74] or they can as well be a portal to Dark matter (DM) [75,76]. Moreover ALPs have been extensively studied in the area of flavour physics, collider physics, astrophysics, cosmology (for extensive discussions see [77,78] and references therein). ALPs being a pseudo-Nambu Goldstone states, their mass can be much smaller than the underlying U (1) global symmetry breaking scale that suppresses their interactions with SM sector, i.e. m a << f a . Thus by taking a large value of f a , one can still make ALP scenario feasible to be explored in multiple current and upcoming experimental facilities (see for example [77]).
Based on the aforementioned discussions, one study the phenomenology of ALPs at the electroweak scale in the effective field theory (EFT) framework. While doing this one is agnostic about the details of the underlying ultra-violet (UV) complete physics at the high scale Λ ∼ 4πf a [79,80]. In this picture, we assume that the ALP a is the SM gauge singlet pure pseudo-scalar with only CP conserving couplings with the SM particles and its mass m a v SM . Based on this ALP EFT framework can play a crucial role in many interesting particle physics phenomenology. For example, ALPs have been suggested to explain the current discrepancy between theoretically expected and measured value of the anomalous magnetic moment of Muon [81],recently claimed apparent resonance in a particular nuclear transition of 8 Be in the Atom ki pair spectrometer experiment [82] and more recently observed anomalous resonance in 4 He [83]. In addition, many collider studies have been carried out to look for the signatures of both stable and decaying ALPs [80,[84][85][86][87]. Motivated by these active research on ALPs, in this paper we plan to explore ALP's role as a portal to FIMP like fermionic dark matter. For this we first introduce a fermionic dark matter χ which is SM gauge singlet in the particle spectrum. The stability of the DM (χ) is obtained by imposing a discrete Z 2 symmetry under which χ is odd and rest of the particles are even. Based on the aforementioned discussed EFT formalism, one can express interaction of the ALP with DM, which makes ALP to act as the portal connecting weakly interacting DM with SM particles. It should be noted that the phenomenology of this ALP-mediated DM [88] is very similar to that of pseudo-scalar portal DM [89,90]. In a nut-shell we in-vestigate ALP-portal to FIMP DM on one hand and identify the parameter spaces that are constrained from current limits from ALP searches and also show the regions which will be probed in upcoming facilities in ALP laboratories as discussed before that also satisfy the correct DM relic via ALP-portal.
In the early universe after the end of inflation, the dynamics of thermal DM production and its decoupling from the thermal bath or non-thermal DM production are highly dependent on the evolution history of our universe. In the standard cosmological framework, one assumes that from the time scale of end of the inflation till the onset of the Big Bang Nucleosynthesis(BBN) at the temperature T BBN ∼ few MeV the universe is mostly radiation dominated (RD) [91,92]. It is interesting to note that there is no underlying guiding principle that can forbid one to assume the possibility of non-standard source that dominates the total energy budget of the universe at early epoch i.e. up to pre-BBN temperature. The presence of such non-standard species may significantly modify the dark matter dynamics at the early universe leading to many interesting consequences on DM phenomenology [93]. It turns out that if the total energy budget of the universe is governed by some non-standard species, the Hubble parameter (H) at any given temperature is always larger than the corresponding value of the Hubble parameter in the standard cosmology at the same temperature.
Thus in this non-standard cosmological picture having larger Hubble parameter, the universe expands faster at earlier time (higher temperature) and eventually at some later epoch (before T BBN ) standard radiation density (ρ rad ) overcomes the non-standard energy density (ρ NS ) of the universe. This phase of the universe is identified by the temperature T eq , where ρ NS (T eq ) = ρ rad (T eq ) [93]. From this equality condition, one can assume that when the temperature T > T eq , the universe was in the non-standard cosmological era. Several studies have shown that WIMP or FIMP relic estimations may differ by several orders of magnitude in non standard cosmological scenario [48,[93][94][95][96][97][98][99][100]. It is observed that freeze-in production of DM is drastically suppressed in a faster expanding universe than in RD universe and hence larger couplings (between the visible and the dark sector) are required in order to obtain observed relic density [48,93,98]. In this paper we allude to this alternative cosmological pathway that may give rise to observed DM relic abundance for larger couplings with the visible sector, and thus boost the detection prospects for freeze-in DM. In this work also we show that assuming a non standard cosmological scenario prior to BBN results a better testability of our ALP portal DM framework in ALP search experiments.
The paper is organized as follows. In Section II we introduce our model and discussed all higher dimensional operators that are relevant for the DM phenomenology. In Section III we discuss the basic formalism for DM freeze-in production via ALP. In Section IV we summarize constraints on ALP-DM parameter space from DM relic density. We also show the effect of reheating temperature (T RH ) on the parameter model parameter space. Section V contains discussions on the DM freeze-in production non standard cosmology, while in Section VI we study the consequences of non-standard cosmology in the numerical estimation of DM relic density. We summarize various laboratory and astrophysical searches of axion like particles in Section VII. In Section VIII we briefly mention the direct detection prospect of DM in this scenario. Finally, we provide our conclusions in Section IX. Some important mathematical results are discussed in the Appendices A and B.

II. BASIC SET-UP
In this current frame work we extend the SM particle content by adding a SM gauge singlet Majorana fermion χ and an axion like particle (ALP) a which is a pseudo Nambu-Goldstone boson of a spontaneously broken global U (1) symmetry at some energy scale Λ which is much higher than the electroweak scale v. We also impose an additional Z 2 symmetry under which χ is odd and rest of the particles are even such that χ is stable and play the role of dark matter candidate in this set up. At low energy, interaction of ALP (a) with SM fields and dark matter (χ) is described by the following Effective Lagrangian up to dimension 5: where, L ALP , L ALP−SM and L ALP−DM denote the effective Lagrangians involving ALP only, ALP with SM fields and ALP with dark matter respectively. The couplings of a to the SM fields and DM arise after integrating out heavy fields at scale Λ ∼ 4πf a producing interactions through dimension 5 operators [79,80]: The Lagrangian containing DM χ can be written as: where, M χ denotes the mass of χ and Cχ fa is the effective coupling of DM and a. The last term in the above equation is crucial for our analysis as it describes the interaction between DM and ALP. The interactions mentioned in eq.(3) along with eq.(4) provides a viable connection between SM and DM via ALP portal, and lead to the following interactions: • Freeze out from SM bath mediated by a e.g. gg ↔ χχ, f f ↔ χχ, f f h ↔ χχ, V V ↔ χχ (where, V ≡ W ± , Z) as well as elastic scattering like f χ ↔ f χ,gχ ↔ gχ, V χ ↔ V χ.
• Freeze out from ALPs when ALPs may or may not be in thermal bath e.g. aa ↔ χχ.
• Freeze in from ALPs when ALPs may or may not be in thermal bath in thermal bath e.g. aa → χχ Detailed description of aforementioned scenarios are discussed in ref. [76,101]. Since ALP-SM and ALP-DM couplings are naturally suppressed it motivates us to study DM freeze-in production from the SM bath with ALP as the mediator. We show the parameter space where both ALPs and χ were not in thermal bath due to their feeble couplings. From the above mentioned Lagrangian in eq.(3) & eq.(4) it is understood that DM can be produced from processes like: ff /gg → χχ, and also from the fusion of other vector bosons. However, the same Lagrangian can lead to different two step production mechanisms like ff /gg → aa followed by aa → χχ. This mechanism of DM production is often named as "sequential freeze-in" in literature [101]. However, we restrict our parameter space where DM relic contributions from the processes like ff /gg → aa with aa → χχ are negligible compared to direct freeze-in (mediated via ALP) 1 . The sequential freeze-in is less effective from the perspective of laboratory searches (for ALP as the mediator) of our scenario which will be addressed in the later part (sec.IV) of this paper. The independent parameters of our scenario are: In principle all five operators may be present at the same time in the DM production process.
Simultaneous presence of several such operators not only complicate the numerical analysis but also hide the crucial underlying dynamics of DM production mechanism we allude to explore. Hence for simplicity we will consider one effective operator at a time describing the ALP-SM interaction to understand the dependencies of the involved parameters better. We will focus in details on the following two cases: 1. Gluon Dominance: in this case we consider DM production only from gluon fusions Here the freeze in production channel will be gg → χχ only.
2. Fermion Dominance: in this case we consider DM production from fermions only Here the production channels will be ff → χχ and ff → χχh. Before electroweak symmetry breaking (EWSB)(T EW ) only ff → χχh process is available as understood from eq.(3). Another production channel will be ff h → χχ which is more phase space suppressed compared to ff → χχh.
The Feynman diagrams are shown in fig.1. DM production from other vector boson fusions are perfectly possible. However, in this work we study only the aforementioned two possibilities and confine our analysis in more minimal set up. Ref. [76] has also studied DM freeze in from SM bath via ALP mediated processes and considered only fermion dominance 1 In the case of direct freeze in from gluon fusion DM relic is proportional to (C G /f a ) 2 .(C χ /f a ) 2 , whereas in the case of "sequential freeze-in", DM relic scales as (C χ /f a ) 4 . Thus "sequential freeze-in" dilutes the connection of DM production and late time observables (which rely on ALP-SM couplings not on the ALP-DM coupling) [76,102]. To avoid such scenario we choose C G /f a , C f /f a C χ /f a and neglect "sequential freeze-in". scenario. In this work we consider gluon dominance scenario also and discuss the prospects of laboratory searches in the non-standard cosmological frame-work which is absent in the existing literature.

III. FREEZE-IN DARK MATTER
From our discussions in the previous section, it is evident that the freeze-in production of DM and DM relic density in the acceptable range can be induced by various dimension 5 operators. This mechanism of freeze in by higher dimensional operators is popularly known as ultra violate (UV) freeze-in [33]. In this case, DM relic density not only depends on various model parameters, but also directly proportional to the reheating temperature (T RH ) of the visible sector [33]. Hence for the two dominant higher dimension operators the DM relic density depends on these model parameters: With these set of parameters in hand, we now discuss the freeze-in production of χ for these two cases one at a time: • Gluon Dominance process • Fermion Dominance process In the first case DM (χ) is produced non-thermally through ALP mediated gluon fusion process. To find the number density of χ produced in this process we solve the following Boltzmann equation,ṅ where, c[F g ] is the collision term corresponding to this process and is given by where, Π i (i = g, χ) denote the Lorentz invariant phase space, with initial and final state 4 momenta p 1 , p 2 and p 3 , p 4 respectively. The dynamics of the process is clubbed in the matrix element |M | 2 gg−χχ and F g stands for the gluon distribution function.As gluons are relativistic particles one should use the relativistic treatment and take F g as the Bose-Einstein (BE) distribution in principle. However,it turns out that at high temperature the thermal averaged cross-section differs by only factor of about 3 if we use Maxwell-Boltzmann (MB) statistics instead of BE and leads to same number at low temperature as shown in Appendix A. Since both BE and MB distribution give more or less same number density, we can use MB distribution for gluon to find the thermal averaged cross-section. Converting the number density (n χ ) to dimensionless variable Y χ (where ,Y χ = n χ /s and s is the entropy density), popularly knwon as the DM yield we get the following equation: where, Y eq g is the equilibrium abundance of gluon and σv gg→χχ is the thermally averaged cross section σ gg→χχ given in eq.(A4). s is the co-moving entropy density and H is the Hubble parameter that decides expansion rate of the universe. The details of yield calculation is given in Appendix A. However using simple dimensional arguments one can see dark matter yield, Y χ has the following functional dependence on various model parameters [33]: In Fig.2 we display the variation of the co-moving abundance Y χ with 1/T for a representative value of C χ /f a = 10 −12 , C G /f a = 10 −10 , M a = 0.1 GeV and M χ = 10 3 GeV with T RH = 10 9 , 10 8 and 10 7 GeV depicted by the red,magenta and green lines respectively.
From Fig.2 we note that with an increase in T RH the freeze-in takes place at earlier time with higher abundance. This behavior can be easily understood from the dependence of Y χ on T RH in eq.(10).
Similar kind of analysis is performed with fermion dominance for which the collision term C[F] in eq.(7) involves processes like ff → χχ and ff → χχh. As realized from eq.(3), that before the electroweak symmetry breaking (EWSB)(T EW ) only ff → χχh process is available, and for T T EW the abundance, Y χ has the following functional form [33,76]: where, N c is the color factor associated with quarks. The detailed numerical formulation is shown in Appendix eq.(B7). Comparing eq.(10) and (11) one notices that for fermion dominance Y χ has same dependence on C f as on C g in gluon dominance apart from the Yukawa coupling and the color factor. We also have considered other production channels like ff h → χχ, f h → f χχ and Y χ will have similar form like in eq.(11). Nevertheless for fermion dominance also, DM will be produced from the UV freeze-in and will have the same functional dependence on T RH as shown in Fig.2. For T RH < T EW , ff → χχ process also becomes relevant. However we find that for low T RH the DM remains under abundant (as Y χ ∝ T RH ) for most of the parameter space. This motivates us to consider T RH > T EW throughout our analysis.
Since we are interested in freeze-in production of DM, the criterion of non-thermal pro-duction of DM χ gives us an upper bound on the interaction rates. Hence, to forbid the thermalization of χ one must satisfy the following conditions, where, H is the Hubble expansion rate (H ∼ 1.66 T 2 /M pl ) and Γ gg→χχ is the interaction rate of the process gg → χχ given by n eq g σ gg→χχ v where n eq g is the equilibrium number density of gluons. Similarly Γ ff →χχh is the interaction rate for ff → χχh given by n eq where n eq f is the equilibrium number density of fermions. In this UV freeze in scenario, since most of the DM production happens at a temperature T = T RH , we check the DM thermalization condition at T = T RH . As the interaction rate increases with increasing temperature, the annihilation process might get thermalized at high T RH . This leads to an upper bound on T RH beyond which χ will be thermalised and the DM mass M χ upper limit depends on the couplings C g (C f ), C χ and the mass of DM M χ . From eq. (12) and (13) one can approximate the upper limit on T RH as, T where, M pl is the Planck mass and the details are shown in appendix A. From eqns. (13)- (14), the upper limits on T RH can be understood. On the other hand, the lower bound on T RH > few MeV comes from Big Bang Neucleosynthesis (BBN) [91,92].

IV. RELIC DENSITY
In the previous section we have discussed that the DM (χ) abundance is generated due to the non-thermal production from SM bath via ALP portal. The abundance of χ is decided by the annihilation cross section of gluons (fermions) or in other words, by SM-ALP coupling (C G /f a , C f /f a ) as well as ALP-dark sector coupling (C χ /f a ) as understood from our earlier discussions. We estimate the relic density of χ after solving for Y χ from eq.(7) and is given by, In this UV-freeze-in process besides the couplings and M χ , the relic density also depends on T RH . We are interested in ALP mass M a ∼ O(MeV-10 GeV) as this is this the ballpark region that is suitable to long-lived ALP-portal searches as we will show in subsequent section. However the freeze-in production takes place mostly at T RH ∼ a few TeV which is much higher than the mediator mass. In this limit the thermally averaged cross-sections, respectively (see eqs. (14) or (15)). In Fig.3 we present the (non-)thermalization constraint in C G /f a vs C χ /f a plane for gluon dominance and C f /f a vs C χ /f a plane for fermion dominance. Within the wide range of allowed parameter space the relative strength between these effective couplings (C G /f a , C f /f a and C χ /f a ) will decide the interplay of SM sector, ALP and dark sector leading to numerous production mechanisms of DM as discussed in Ref. [76,101]. However, we restrict our analysis in more simplified set up where we study only DM production from SM bath mediated by ALPs. ALPs produced from SM particles may get thermalized in early universe leading to ALP freezeout from thermal bath making the scenario more complicated. For this reason we restrict our analysis to the parameter space where ALPs do not get thermalized. This helps us to understand the underlying dynamics of the ALP portal dark matter freeze-in process and enable us tp probe the corresponding parameter space after taking into account constraints from various laboratory searches of ALPs. To prevent the production of ALPs we impose the thermalization condition on the processes ij ↔ aa and i ↔ aa (i, j = SM particles) [76].
We incorporate these conditions and the viable parameter space is shown in Fig.3. In the top panel in Fig.3(a) and Fig.3 We discuss the laboratory search prospects of our scenario in sec.VII. However due to very tiny couplings most of the parameter space satisfying the observed relic is not within the reach of ALP search experiments as shown in Ref. [76]. This motivates us to explore the same scenario in modified cosmological framework which will be discussed through out the latter part of this paper.

V. NON STANDARD COSMOLOGY
In this work hitherto we assume the universe at the time of freeze-in production was radiation dominated (RD) whose energy density (ρ rad ) red-shifts with scale factor a as ρ rad ∼ a −4 . The ALP portal couplings between visible and dark sector that provides the correct observed relic density is too small to be explored in current laboratory search for ALPs, as shown in [76,103]. On the other hand the effective coupling of DM pair with ALP also requires to be small to keep the DM and ALP non-thermalized in the early universe (eq.(12) & eq. (13)). . This small ALP-portal coupling automatically leads to highly suppressed DM nucleon scattering cross-section, thus evading current dark matter direct detection limits, which is the general feature of any freeze-in dark matter production process.
The dark matter relic density highly sensitive to the cosmic evolution history. Any drastic change in the dynamics of the cosmic evolution may have nontrivial impact on the dark matter relic density. So far, while discussing the dark matter relic density we have considered that the universe was radiation dominated (RD), i.e. within the standard cosmological paradigm. And this is motivated from the fact that the universe was RD at the time of BBN only [91,92]. However, due to the lack of our understanding of the nature of universe before BBN, one can possibly consider that our universe was not radiation dominated prior to BBN. Thus it is fair to assume that at early epoch of time, the universe was dominated by some non-standard species and the corresponding scenario is popularly known as the non-standard (NS) cosmological framework [94,98,104]. The DM abundance in freeze-in scenario changes rather dramatically if we assume cosmic history other than RD at the time of freeze-in [94,98]. realized shortly at the end of this section. We present a brief summary of NS cosmology here. Non standard cosmology may be realized by introducing a scalar field φ with equation of state (EOS) p φ = ωρ φ , where p φ , ρ φ denote the pressure and energy density of φ respectively and ω is the EOS parameter [104]. The energy density of such a species dilutes with scale factor z as where, n is given by n = 3ω − 1. For n > 0 i.e. ω > 1 3 , ρ φ red-shifts faster than radiation energy density [105][106][107][108][109].
In the early stage or at high temperature, in presence of this φ the total energy density (ρ tot ) of the universe becomes ρ tot = ρ rad + ρ φ . As the universe cools down, it becomes radiation dominated again prior to BBN temperature (T BBN 1 MeV) [91,92], so that the BBN predictions of light element abundances remain unaltered and during this process at some temperature T eq T BBN , ρ rad becomes equal to ρ φ . If one considers T eq close to T BBN 1 MeV, the number of active neutrino flavours (N ν ) receives an additional contribution from ρ φ which can significantly modify bound on ∆N ν ≡ N ν − N SM ν from the BBN [110]. Thus to comply with this BBN bound T eq and n are restricted by the following condition [104]: Following ref. [104] using the definition of T eq and entropy conservation one can write the total energy density as, where, g * and g * s are the relativistic degrees of freedom associated with energy density and entropy density. With this total energy density in NS cosmology the Hubble parameter is defined as H(T ) ∼ √ ρ tot /M pl in contrast to RD universe where Hubble parameter scales as H(T ) ∼ √ ρ rad /M pl . Since for any n > 0 and T > T eq , ρ tot is larger than ρ rad ,H(T ) is also larger in NS cosmology compared to the case when the universe is only radiation dominated.
Thus NS cosmology leads to the scenario when the universe expands faster than the usual RD universe. At the limit T T eq the Hubble scales as [104] We will now study the phenomenological consequence of non standard cosmology in the context of our present set up and show how the analysis differ from that in RD universe for any n > 0.
After introducing the basic idea and motivation of non standard cosmology we are now set to discuss the modified freeze-in abundance of DM in NS cosmological era. In such fast expanding universe, to keep the DM particle χ non-thermalized at early universe, eq. (12) has to be satisfied but with a modified Hubble parameter as in eq. (20). In Fig.4 (19). From this figure it is also evident that with higher values C G /f a ( 10 −3 ), the process gg → χχ gets thermalized in RD universe (n = 0) but fails to thermalize in non-standard (NS) scenario (n = 0). This leads to the relaxation of the parameter space in the effective coupling plane (C g /f a vs. C χ /f a ), allowing larger values of these effective coupling for the non-thermal production of DM which will be analyzed in the next section. Similar kind of conclusion can be also drawn for fermion dominance also i.e. for the evolution of Γf f →χχ .
To understand the impact of this faster than radiation expansion history of the universe on DM abundance we again solve eq.(7) with the Hubble parameter, H defined as in eq. (20).
We now have two more free parameters n and T eq in addition to what we had in standard cosmological scenario. The free parameters in NS scenario are given by, Considering gluon dominance we solve eq.(7) to evaluate the co-moving abundance Y with and depicted by black, magenta, red, blue, green respectively. From the above mentioned plot we notice that for a fixed T eq higher value of n results lower value of Y . From eq. we perceive that for same n > 0 with higher T eq leads to higher DM abundance. This can be readily comprehended by noting the functional dependence of H on T eq in eq.(20) which apprises that an increase in T eq leads to a decrease in H consequently an increase in Y . In the aforementioned figure, the line corresponding to n = 0(black line) which represents RD universe remains unaltered with the variation of T eq . This is also true for fermion dominance case.
From the above discussion we conclude that same couplings will produce lower abundance in fast expanding universe than in RD universe. Needless to say while other parameters are kept fixed, one has to pick higher values of effective couplings in fast expanding scenario than in RD universe in order to satisfy correct relic density without violating the thermalization bound. Such large value of effective ALP couplings can be further probed using various laboratory searches of ALP. This motivates us to explore the ALP phenomenology in the non standard cosmological scenario in the rest of this paper.

VI. RELIC DENSITY IN NON-STANDARD COSMOLOGY
Following the discussion of the previous section here we investigate the consequences of non-standard cosmology in deciding the relic density. In Fig.6 we present the allowed parameter space where the DM production processes do not thermalize at the time of freezein. We show the parameter space in C G /f a vs C χ /f a plane for gluon dominance and C f /f a vs C χ /f a plane for fermion dominance. In Fig.6(a) and Fig.6 From all four plots in Fig.6 we notice that with increase in n, the allowed parameter space also increases in the effective coupling plane. This feature can be explained using eq. (20) where we see the expansion rate of universe H increases with increase in n and with this larger value of H the process fails to thermalize even with higher values of C G /f a (or, C f /f a )and C χ /f a satisfying eq. (12). Comparing plots in left panel ( Fig.6(a),6(c)) and right panel ( Fig.6(b),6(d)) one notices that with lower value of M χ the allowed parameter space is larger. The reason behind this is that the interaction rate is proportional to M 2 χ (eq. (A8)) and smaller M χ results smaller interaction rate. Comparing Fig.3 and Fig.6 one notices that for a same T RH non standard cosmology provides larger values of effective couplings within as we have shown in the above mentioned figure.
In Fig.7 we present the relic satisfied contours in T eq vs. n plane to showcase the dependence of n and T eq in deciding relic density while other parameters are kept fixed. Here we consider gluon dominance and take fixed values of T RH = 10 7 GeV, C G /f a ∈ (4 × 10 −9 , 4 × 10 −8 , 4 × 10 −6 , 4 × 10 −4 ) shown in magenta, blue, red and black curves respectively. The gray shaded region is excluded from the bound on T eq from BBN as mentioned in eq. (18). One can see that for a fixed T eq higher C G /f a demands higher values of n to satisfy the relic density as we have already pointed out in the context of Fig.5. Hence we see C G /f a = 4 × 10 −4 (black line) contour lies on the top while C G /f a = 4 × 10 −9 (magenta line) in the bottom in T eq − n plane. Another noteworthy feature from the above mentioned figure is that for a fixed value of C g /f a an increase in T eq calls for an increase in n also. This can be explained by the fact that higher T eq results lower value of H and to satisfy the relic density n has to be increased to compensate the effect of T eq .

VII. LABORATORY AND ASTROPHYSICAL SEARCHES
In this section we will discuss the detection prospects of ALP for our set up. We present all sources of astrophysical, cosmological and collider constraints depending on the mass of ALP and effective ALP-SM couplings. Our results show that significant portion of the parameter space satisfying observed relic density in NS cosmology is excluded by various existing constraints. However some of the parameter space satisfying observed relic density in RD universe is far below the sensitivity reaches for all kinds of constraints. We now discuss implications of those constraints on our scenario.
• Astrophysical constraints: Among the astrophysical bounds SN1987A puts the most stringent bound on the effective ALP-SM coupling. In standard picture of corecollapse super nova, it cools by emitting most of its energy via neutrino emission.
However it might get affected in the presence of additional weakly interacting particles which can carry away additional energy from the core. Likewise in our set up ALPgluon couplings (C G /f a ) or, ALP-fermion couplings (C f /f a ) will be responsible for additional cooling via nucleon (N ) scattering (N N → N N a) [111][112][113]. For a large ALP couplings although ALP would be produced in a large number due to stronger interaction, the same coupling will lead to re-absorption of ALPs in the core [111]. For small ALP mass (M a < 100 MeV) the constraint from SN1987A is relevant and we incorporate that bound in our analysis [103,111]. There is another relevant constraint coming from Horizontal Branch (HB) star [114,115]. Any pseudo-scalar that can couple to electron or photon can interact with particles in the core of HB stars and disrupt the stellar evolution. Similar to SN1987A , ALPs can also affect the cooling mechanism of HB stars and hence constraints from HB stars put limit on the effective coupling of ALPs [76]. However this bound is sensitive only in the keV range of M a [114].  [118]. Here ALP can be produced by bremsstrahlung or positron annihilation and decay to visible particles (a → e + e − , γγ) via either tree level or loop level processes. On the other hand for small effective coupling with SM, ALP becomes a long lived particle (LLP) which can be separable from background [119]. These ALP's can be produced in proton collision and can be within the reach of displaced decay search with LHC Track trigger [103]. The FASER experiment at LHC will also be able to probe ALP's with mass above few MeV [120]. Another bound on ALPs come from the decay of heavy mesons and this bound is sensitive in M a ∼ MeV to few GeV region. For processes like B + → K + a(→ µ + µ − /e + e − ) LHCb puts limit on effective ALP-SM coupling [121]. Experiment like CHARM also constrains the effective coupling considering the process B 0 → K * 0 a(→ µ + µ − /e + e − ) [122]. For MeV scale mass another constraint arises from NA62 which is designed to measure rare Kaon decays. NA62 puts constraint on the effective ALP-SM coupling using the decay of charged Kaon(K + ) to pion(π + ) and ALP, where ALP acts as the long lived particle and [111]. However this experiment is sensitive to M a < M π and loses sensitivity for M a ≈ M π due to large background. DUNE (Deep Underground Neutrino Experiment) can also be used to search for ALPs where ALPs are produced with meson mixing or with gluon fusion [84] and later they decay to photons or hadrons in the near detector (ND) [103].
After demonstrating all kinds of existing constraints and future sensitvity reaches we present the result in M a vs. effective coupling plane in the following two subsections.

A. Gluon dominance scenario
Here, we consider gluon dominance scenario where only C G /f a and C χ /f a non zero while rest of the effective ALP couplings are zero. We present all the constraints discussed above in the M a vs. C G /f a plane. We also portray the parameter space that reproduce observed relic density in the same plane.
In Fig.8 we showcase all the constraints as well as the future experimental reaches on  Fig.8(a) and T eq = 2 GeV in Fig.8 shown by blue, green and red dashed line respectively in the top panel. From these two plots we notice that the higher the value of n the better the experimental sensitivity as discussed in sec.VI. We also note that for RD (n = 0) universe and for n = 1 (NS) the relic satisfying contours (C G /f a < 10 −10 ) are far below the existing and future constraints and we do not portray them here. This is what we expect from the thermalization condition discussed in sec-IV. Comparing Fig.8(a) and Fig.8(b) we notice that for a fixed n, the lower the value of T eq the higher the value of C G /f a required (as discussed in sec.VI) and hence better experimental sensitivity. From the above figures we notice that though in RD universe the relic satisfying couplings for M χ = 10 3 GeV are far below experimental sensitivity but significant amount of parameter space reproducing observed relic density is excluded in NS cosmology. SN1987A and BBN puts most stringent bound for M a < 100 MeV and lower values of C G /f a . The region "Existing Constraints" include the bounds from partially invisible kaon decays [111], electron beam dump, CHARM, visible kaon decays, B decays [123], LHC dijet searches [124]. However, there is still some unexcluded parameter space in NS cosmology which can be probed in future experiments like DUNE ND (C G /f a > 10 −9 ), LHC track Trigger (C G /f a > 10 −8 ) and FASER(C G /f a > 10 −6 ).
In Fig.9 we showcase all the constraints as well as the future experimental reaches in M a vs. C G /f a plane with M χ = 1 GeV and other parameters are kept same as in Fig.8. The horizontal colored lines in black,magenta, blue dashed line the contours satisfying observed relic density for different values of of n ∈ {0, 1, 2} with T eq = 20 MeV in Fig.9(a) and T eq = 2 GeV in Fig.9(b) respectively. From these two plots we again opine that the higher the value of n better the experimental sensitivity as discussed in sec.VI. Relic satisfying lines for n ≥ 2 in Fig.9(a) and n ≥ 3 in Fig.9  in Fig.10. The horizontal colored lines in Fig.10(a) & Fig.10(b) are the contours satisfying observed relic density with fixed values of T RH = 10 6 GeV and C χ /f a = 10 −10 . We showcase the relic density contours with different colored lines for different values of n with T eq = 20 MeV and T eq = 2 GeV in Fig.10(a) & Fig.10(b) respectively with fixed value of M χ = 10 3 GeV. The contours with different values of n ∈ {0, 1, 2, 3, 4} are shown by black, magenta, blue, green and red dashed line respectively. From these two plots we observe that higher values of n give better experimental sensitivity keeping other parameters fixed, as discussed in sec.VI. Here also we claim that for RD (n = 0) universe the relic satisfying contours (C f /f a < 10 −10 ) are far below the existing constraints as expected from the thermalization condition discussed in sec-IV. From this plot we conclude that for a fixed n lower values of T eq provides better experimental sensitivity than higher values of T eq (as discussed in sec.VI).
Similar sets of constraints are obtained in Fig.11 Inspite of all of these stringent limits, there is still some unexcluded parameter space providing observed relic density in NS cosmology which can be probed in future generation experiments. In Fig.12 we present such future projection along with relic satisfying contours in M a vs. C f /f a plane. There is a mass window for M a = 0.2 − 1.3 GeV which can be explored in future as discussed in detail in ref. [125] and we present the same in the context of our analysis. Here in Fig.12 also the horizontal colored lines in all plots depict contours satisfying observed relic density with fixed values of T RH = 10 6 GeV , C χ /f a = 10 −10 and M χ = 10 3 GeV. We exhibit the relic density contours with different colored lines for different values of n with T eq = 20 MeV in Fig.12(a) and T eq = 2 GeV in Fig.12(b) respectively with T eq =20 MeV LHCb: B→ K * a K + → π + a(→ µ + µ -) K + → π + a(→ e + e -) E137(→ e + e -) E137 (→ γ γ) K + →π + inv HB star SN1987A

T eq =2 GeV
LHCb: B→ K * a K + → π + a(→ µ + µ -) K + → π + a(→ e + e -) E137(→ e + e -) E137 (→ γ γ) K + →π + inv HB star SN1987A in Fig.13(a) and T eq = 2 GeV in Fig.13(b) respectively, while other parameters are taken same as in Fig.12(a) and Fig.12(b) respectively. For M χ = 1 GeV in Fig.13(a) and 13(b) we don't show the relic satisfying lines for n ≥ 2 as the required coupling C f /f a > 10 −3 is already excluded from existing constraints. For brevity we don't discuss again the features of the contours satisfying relic density here. However, it is easy to infer that, while holding all other parameters constant, contours characterized by higher values of n and lower values of T eq are more sensitive to future projections. Likewise when n and T eq are held fixed, we find that sensitivities of future experiments are such that lighter DM masses can be easily probed.
The gray shaded region in all the plots in Fig.12 and Fig.13 is excluded from ArgoNeuT experiment at Fermilab where they serach for ALPs produced from neutrino beam target and deacying to dimuon pairs [126]. Experiments like SBND [127], ICARUS [127], SHiP [128] and FASER 2 [129] will be able to probe a significant amount of parameter space reproducing observed relic density in NS cosmology. The label "ICARUS * " signify ICARUS experiment with off axis neutrino beam at the Main Injector [125]. The label "DUNE+" depicts the result of DUNE experiment with additional 30m length before the front detector panel as part of the decay volume as discussed in ref. [125].
Therefore we demonstrate that our proposed frame work is more testable in future gen-eration experiments 3 for both gluon dominance and fermion dominance scenario when the fast expansion parameter n takes higher values. This was the primary motivation behind this work.

VIII. DIRECT DETECTION
To get signal in direct detection the DM in our set up χ has to scatter nucleon. The elastic scattering can be mediated by the effective Lagrangian where, we assume the momentum transfer q is much smaller than M a . From the effective interaction in the above equation we notice the presence of γ 5 in the χ χ a vertex which will lead to a suppression of q in non-relativistic limit in the amplitude level [131]. The a G G vertex will also reduce to a current containing γ 5 which will eventually lead to suppression of another factor q [132,133]. So the overall cross-section will contain a suppression of q 4 which is a general feature of direct detection mediated by pseudo-scalar [134][135][136]. Apart from the small effective couplings (C G /f a , C χ /f a ) required to produce χ non thermally the direct detection cross-section gets a suppression of q 4 where as typical value of q 2 in current experiments are O(10KeV) [16] and hence evade the current experimental sensitivity.

IX. DISCUSSION AND CONCLUSIONS
In the parameter space for ALP decay constant not satisfying the limits f a ≤ 10 9 − 10 12 GeV, the axions cannot behave as the dominant cold DM component of the universe via the traditional misalignment mechanism. However ALP may still act as portal to DM; in this work we have studied a scenario where an ALP act as a portal to dark sector in an effective theory framework in radiation-dominated early universe and other non-standard cosmological histories and proposed detectable regions in such parameter space to be within the reach of future ALP search experiments satisfying the observed DM relic density. From our investigations we summarize our major findings below: • DM relic density is proportional to the effective couplings C G /f a ,C f /f a ,C χ /f a and M χ (mass of DM) as shown in eq.(10) and eq. (11). In our analysis we consider higherdimensional non-renormalizable operators between dark sector and visible sector via the ALP-portal, therefore the DM production have dependencies on T RH (see eq.(10), (11)). We find that the relic density is also proportional to T RH (see Fig.2). However, for our chosen parameter space the relic density is almost insensitive to the light mediator M a due to freeze in happening at high T RH . Conditions from the nonthermalisation of χ sector with the SM bath imposes a strict limit from the DM production on the effective couplings as decided by M χ dependencies in eq.(14) and eq.(15). Our analysis shows that for M χ 50 MeV the effective couplings are already excluded, whereas M χ ∼ 1 TeV the effective couplings are too small to be probed by the present and future generation experiments.
• The aforementioned situation changes drastically if we assume a non standard species φ contributing significantly to the energy density of early universe at an early epoch. In such scenario the expansion rate of the universe H is modified which leads to a faster expansion of the universe than the radiation-dominated universe. In such a framework the value of H allows larger values of effective couplings (C G /f a ,C f /f a ,C χ /f a ) still keeping the dark sector non thermalized. The relaxation of parameter space can be understood by comparing Fig.3 and Fig.6 where we show that the effective couplings can be O(10 2 ) larger in NS cosmological scenario than in RD universe. After solving respective Boltzmann equations with modified Hubble H we are able to identify the parameter space satisfying the observed relic density. This boosts the detection prospects for the scenario in laboratories.
• The non-standard cosmology allows the involved effective SM-ALP couplings to be large in order to satisfy the observed relic, therefore this can be probed in present and future generation experiments. Moreoever from Fig.8 and Fig.10 we find that higher the value of n, the better the detection prospects for ALP-SM couplings. The most stringent bound on the relic satisfying parameter comes from astrophysical probes specially SN1987A and HB stars. The value of ∆N ef f (dark radiation) estimated from BBN gives a lower bound on the ALP mass M a . Inspite of taking into consideration the existing bounds on C G /f a and C f /f a coming from electron beam dump, CHARM, visible kaon decays, B decays and LHCb, etc. There remains still leaves some unexcluded parameter space satisfying the observed DM relic density. The DUNE and LHC track trigger will be able to probe such experimental parameter space that coorespond to C G /f a ∼ 10 −5 − 10 −9 for the mass range of M a ∼ 0.1 − 10 GeV as shown in Fig.8. Similarly, experiments like SHiP, DUNE, FASER2, ICARUS will be able to probe C f /f a ∼ 10 −4 − 10 −6 for the mass range of M a ∼ 0.2 − 1 GeV as displayed in Fig.12.
• Finally we point out that our analysis is also an alternative probe for non-standard cosmological evolution of early universe. In Fig.7 we show that for a fixed (C G /f a ), (C f /f a ) and C χ /f a one is probing values of n and T eq which satisfy the observed DM relic density. The detection prospects of various effective couplings for our proposed scenario are also highly sensitive to n and T eq . For example, comparing Fig.8(a) and 8(b) we note that the future probable regions of relic satisfying contours change with T eq for the case of gluon dominance. Similar conclusion can be drawn for the fermion dominance case by comparing Fig.12(a) and 12(b). Hence any signal in the upcoming experiments will provide us concrete information regarding the cosmic evolution and dark sector assuming ALP-portal freeze-in DM scenario.
Although we have several possible ways by which the parameter space for ALP-mediated For the case of Gluon dominance scenario, freeze in production of DM from SM bath takes place via the following two operators : Including the color factor the amplitude square will be, |M | 2 gg−χχ = 256 where, s is the center of mass energy. The cross-section will be Consider the process: g + g → χ + χ The Boltzmann equation is given by, n χ + 3Hn χ = dΠ g dΠ g dΠ χ dΠ χ δ 4 (p 1 + p 2 − p 3 − p 4 )|M | 2 gg→χχ F g (p 1 ) F g (p 2 ), where, p 1,2 are the momentum of incoming gluons and p 3,4 are the momentum of outgoing χ particles. F g (p i ) denotes the distribution function of gluons.

Using MB distribution:
We assume initial bath particles follow the MB distribution, F g ∼ e −Eg/T . The term in the R.H.S of eq (A5) is the collision term C[F g ]. After some simplification the collision term looks like [33]: where, P ij is given by, Upon proper substitution in eq. (A6) we get,

Using BE distribution
Using the proper quantum correction the collison term takes the following form [138], where, F g is now BE distribution function and we write it as [138], where, V µ is the 4-velocity in the center of mass frame The cross section for the process gg → χχ is given by, where E 1,2 initial sate particle's energy and v mol is the Moller velocity. v mol is given by, So, incorporating v mol in eq.(A11) one obtains the following equality: Now we define the collison term as C rel [F g ] = dΠ 1 dΠ 2 4p 1 p 2 σ gg→χχ v mol F g (p 1 )F g (p 2 ) Following the notation in ref. [138], in gas rest frame we define where, p is related to c.o.m frame by a Lorentz transformation p = Λ(E, 0, 0, 0) , E being the c.o.m energy and Λ is the Lorentz transformation. Λ is function of rapidity η and angular variables θ, φ as given in [138]. Similarly k can also be translated to c.o.m frame k = Λ(η, θ k , φ k )(0, l 1 , l 2 , l 3 ), where |l i | = | p 1 − p 2 |. where θ k , φ k are angular variable associated with k-space. So, for mass less particle |l i | = E. For convenience we consider the momentum along only one direction as k = Λ(η, θ k , φ k )(0, E, 0, 0).
Following the formalism as in ref. [138] we can write the distribution functions as After little algebra we can also write the collision term in eq.(A14) accordingly [138], where dΩ i ≡ sin θ i dθ i dφ i is the solid angle. We perform the solid angle integration first and get, Incorporating this in eq. (A17) we will find Γ gg→χχ = 1/n eq g C rel [F g ] using BE distribution. In Fig.14 we compare the two interaction rates for MB distribution and BE distribution.
We denote Γ BE and Γ MB as interaction rates for g +g → χ+χ using BE and MB distribution function respectively. We show the evolution of the ratio Γ BE /Γ MB with 1/T (GeV −1 ) for C χ /f a = 10 −12 , C G /f a = 10 −10 using two different values of M χ = 10 3 , 10 4 GeV depicted by the red and blue solid lines. We notice that at higher temperature only, the ratio differs from 1 i.e. the two interaction rates Γ BE and Γ MB differ by a factor of ∼ 3. However, at low temperature the ratio becomes 1. This justifies our approximation using MB distribution.