Thermally corrected masses and freeze-in dark matter: a case study

If coupled \emph{feebly} to the Standard Model bath, a dark matter can evade the severe constraints from the direct search experiments. At the same time, such interactions help produce dark matter via the freeze-in mechanism. The freeze-in scenario becomes more interesting if one also includes the thermal masses of the different particles involved in the dark matter phenomenology. Incorporating such thermal corrections opens up the possibility of dark matter production via forbidden channels that remain kinematically disallowed in the standard freeze-in setup. Motivated by this, we investigate such freeze-in production of the dark matter in a minimally extended $U(1)_{L_\mu-L_\tau}$ framework that remains consistent with the recent muon $(g-2)$ data. Here, the role of the dark matter is played by the scalar with a non-trivial charge under the additional symmetry $U(1)_{L_\mu-L_\tau}$. This scalar dark matter obtains a thermally corrected mass at high temperatures for a not-so-small self-coupling. We show that the thermal correction to the dark matter mass plays a significant role in the dark matter phenomenology.


I. INTRODUCTION
The existence of a non-luminous and a non-baryonic form of matter in the universe popularly known as dark matter (DM) [1,2] has attracted the attention of the scientific community in the past several decades. Several vital pieces of evidence, like galactic rotation curves, gravitational lensing, and the anisotropies in the cosmic microwave background, have confirmed the presence of this form of matter in the universe. The vital information known about the DM is its relic density which is very accurately measured by the experiments studying the anisotropies in the cosmic microwave background radiation [3,4]. Despite this vital information, the particle nature of the DM is still unknown. On the other hand, the Standard Model (SM) also fails to provide a particle that can be identified as a DM candidate. This suggests that physics beyond the SM is inevitable.
Depending on the interactions of the DM with the SM bath, several theories have been proposed in the literature that tries to explain the particle nature of the DM. Among these theories, the most popular DM candidate is the weakly interacting massive particle (WIMP) [5][6][7][8][9][10][11][12][13][14][15]. Such a DM candidate thermalizes in the early Universe if the temperature of the thermal bath is above its mass. As the Universe cools down and the temperature of the plasma reaches below the mass of the DM mass, its abundance freezes out. Due to its not-so-small interactions with the SM particles, the WIMP type of dark matter is subjected to various experimental constraints. The experiments like LUX [16,17], PANDA [18], XENON1T [19], LUX-ZEPLIN (LZ) [20] provide a stringent bound on the DM-nucleon scattering cross-sections. Such severe constraints can be easily evaded if the DM candidate is a feebly interacting massive particle or FIMP [21][22][23][24][25][26][27][28][29][30]. As the name itself suggest, this kind of DM interacts very feebly with the SM particles and hence cannot be tested in the experiments mentioned above. Due to the feeble interaction, a FIMP-type DM is slowly produced from the thermal bath, and once the number density of the bath particle that is responsible for its production becomes Boltzmann suppressed, the DM abundance freezes in.
Recently, studies like [24,29] have shown incorporating the thermal masses [31][32][33][34] of the species involved in the DM production can open up the kinematically forbidden channels that are not allowed in the standard freeze-in (SFI) scenarios. Here, it is assumed that the particle responsible for the production of the DM is not only a part of thermal plasma but also acquires a significant thermal mass at high temperatures due to its interactions with the thermal bath. As a result of the thermal corrections, the thermal mass of the species present in the plasma can be significantly different from their bare masses. Consequently, even if a DM production channel is kinematically disallowed (forbidden) at T = 0 due to the hierarchy of masses between the bath particle and the dark matter state, that can change at a higher bath temperature. A T = 0 correction to the mass of the bath particle can flip the hierarchy, thereby kinematically opening up that particular production mode.
In the specific context of DM production through the decay of a bath particle that picks up a thermal mass M (T ), the hierarchy M (T ) > 2M DM is obtainable through non-zero temperature effects. Such freeze-in production of the DM via forbidden channels is also known as forbidden freeze-in (FFI). An important study is [35] that correlates cosmological phase transition, a phenomenon also driven by T = 0, to freeze-in in the presence of thermal corrections. Another one discussing thermal effects in DM and gravitational wave signatures is [36].
Unlike the standard U (1) B−L model [11,29,43], which offers a stable DM in the form of right-handed neutrino (RHN), the U (1) Lµ−Lτ scenario requires an additional scalar (singlet under the SM gauge symmetry) with a non-trivial U (1) Lµ−Lτ charge to explain the presence of the DM in the Universe. Assuming the DM interacts feebly with the bath particles, it can be produced from the decay of (i) the scalar responsible for the breaking of U (1) Lµ−Lτ symmetry, (ii) SM Higgs and, (iii) the massive gauge boson of U (1) Lµ−Lτ symmetry in the SFI scenario 1 if kinematically allowed, as was also discussed in [39]. In this work, we aim to explore the deviation observed from the SFI scenario once the thermally corrected masses for the bath particles are taken into account. Besides explaining the dark matter, an U (1) Lµ−Lτ framework can simultaneously explain the discrepancy in the anomalous magnetic moment of muon (g − 2) from its SM prediction [40] and non-zero neutrino masses [40]. Keeping this in mind, we show that the present setup can accommodate a DM that can be produced via both SFI and FFI channels while also providing the solution for the discrepancy in the muon anomalous magnetic moment.
The paper is organized as follows. The model is introduced in section II, and the various constraints deemed relevant are detailed in section III. We compute the relevant thermal masses in section IV. And the same section also elaborates on the ensuing freeze-in phenomenology. Finally, the study is concluded in section V.

II. THE MODEL
We extend the SM gauge symmetry by an U (1) Lµ−Lτ symmetry where L µ and L τ represent the muon and tau lepton numbers, respectively. The fermionic content of the model includes the SM leptons and quarks together with three additional right-handed neutrinos (N e , N µ , N τ ). As suggested by the symmetry of the present scenario, the muon and tau carry a non-trivial charge under the U (1) Lµ−Lτ . The newly introduced RHNs are singlets under the SM gauge symmetry, while two of them carry 1 and −1 unit of U (1) Lµ−Lτ , the third remains uncharged. We remind the readers that the addition of RHNs carrying nonzero U (1) Lµ−Lτ charges is necessary to cancel anomalies. The scalar sector of the setup is enhanced with a complex scalar (S) which is a singlet under the SM gauge symmetry but carries 1 unit of U (1) Lµ−Lτ charge. We also introduce an additional scalar (φ), a SM gauge singlet that plays the role of the DM. The stability of the DM is guaranteed by its non-trivial charge assignment under U (1) Lµ−Lτ symmetry. The gauge quantum numbers of the relevant fields are shown in Table I and Table II.  With an idea of particle content and their charges under the different symmetry groups, we now proceed to write their interactions. To begin with, we first write the kinetic terms for the additional fields, where D δ = ∂ δ + ig µτ Q µτ (Z µτ ) δ with Q µτ representing the charge and Z µτ being the gauge boson of U (1) Lµ−Lτ symmetry. Finally F αβ µτ = ∂ α Z β µτ − ∂ β Z α µτ . We also consider the pure U (1) Lµ−Lτ model where tree-level mixing of Z µτ with the SM gauge bosons is absent. from the µand τ -families. Next, we write the Lagrangian involving the Yukawa interactions and masses of the additional fermions involved, Finally, we write the most general scalar potential involving all the scalars in the present setup, The scalar S breaks the U (1) Lµ−Lτ symmetry once its CP even component develops a non-zero vacuum expectation value (vev) v µτ . The U (1) Lµ−Lτ gauge boson consequently obtains a non-zero mass, m Zµτ = g µτ v µτ . The same breaking also results in an additional non-zero mixing that develops between the RHNs, as seen from Eq. 2. After the Electroweak Symmetry Breaking (EWSB), the Higgs doublet (H) also develops a non-zero vev v = 246 GeV. The scalars after the breaking of the gauge symmetry can be parameterized as, Subsequent to the EWSB, a non-zero h − s mixing leads to the following mass terms, The mass matrix is diagonalized using The mass eigenstates (h 1 , h 2 ) then have masses Finally, after both the symmetries are broken, the dark matter mass can be expressed as, It is convenient to describe a framework in terms of vevs, physical masses, and mixing angles.

III. CONSTRAINTS AND ADDITIONAL ISSUES
We discuss in this section the various constraints on this model, as well as the predictions of neutrino mass and the muon anomalous magnetic moment in this model.

A. Neutrino scattering experiments
Gauging the U (1) Lµ−Lτ symmetry leads to severe constraints from neutrino trident pro- Here N denotes a heavy nucleus.
Given the good agreement of the observed results with the SM for this process reported by CHARM-II [45] and CCFR [46,47], the parameter space in the presence of a new neutral gauge boson gets seriously restricted. A relatively new probe is the coherent elastic neutrinonucleus scattering (CEνNS) process that, at the amplitude level, looks like νN → νN .
CEνNS has recently been measured by the COHERENT collaboration [48][49][50]. The Z µτ in the gauged U (1) Lµ−Lτ model enters CEνNS through kinetic mixing with the SM gauge bosons and thereby gets constrained. In fact, the BOREXINO [51,52] process studies the same scattering process as COHERENT, however, using solar neutrinos. The impact of all the constraints is conveniently depicted in the m Zµτ − g µτ plane in Fig. 1.

B. LHC constraints
Z → 4µ searches by ATLAS [53] and CMS [54] constrains the gauge sector of the model and rules out a portion of the m Zµτ − g µτ plane as shown in Fig. 1. On another side, an h − s mixing as defined by Eq.(6) implies that the tree-level couplings of h 1 with the SM fermions and gauge bosons scale by a factor of c θ w.r.t. the corresponding SM ones. This subjects the mixing angle θ to Higgs signal strength constraints and other exclusion limits from the LHC [55,56]. We adopt |s θ | < 0.1 in this work to comply with all such constraints.
Finally, the reported upper limit on the invisible branching ratio of the 125 GeV Higgs puts a limit on BR(h → φφ † ) whenever kinematically allowed. However, this constraint is almost trivially satisfied in this setup given a feeble h − φ − φ † interaction strength dictated by the freeze-in dynamics.

C. Dark matter constraints
We demand that the freeze-in relic density predicted by this model must entirely account for the observed relic of the universe. The Planck experiment [3] has reported In addition, any dark matter model must abide by the limits imposed by the various direct detection experiments mentioned before and the most stringent constraint currently comes from LZ. Direct detection scatterings in this model are mediated by gauge bosons through kinetic mixing as well as through the scalars h 1 , h 2 . It is, however, much easier to evade such constraints in a freeze-in framework such as the present scenario. This is because demanding the observed relic density through the freeze-in dynamics would necssitate that couplings of the DM to h 1 , h 2 are tiny and the charge n µτ is accordingly small. The direct detection cross sections are thus expected to be below the stipulated bound.
The dominant contribution to ∆a µ in this model comes from the 1-loop amplitude mediated by the Z µτ gauge-boson 2 . This contribution can be expressed as [57,58] ∆a Zµτ where r = (m Zµτ /m µ ) 2 . Following the announcement of the FNAL [59] results on muon g − 2, a combined measurement of the discrepancy is An inspection of Fig The cyan band is the region compatible with the 2σ limit of muon g − 2 as quoted in Eq.13.

E. Neutrino mass
Generation of neutrino mass in the gauged U (1) Lµ−Lτ model occurs through Type-I seesaw [60][61][62][63][64] and has been discussed in detail in [14,15,40]. The light neutrino mass matrix has the familiar Type-I form Here, M D and M R refer to the Dirac and Majorana mass matrices. Following the spontaneous breaking of the gauge symmetry of the model, one derives We refer the reader to [14,15,40] for details of fitting the neutrino data. A similar approach is adopted for this study.

IV. FREEZE-IN PRODUCTION AND THE IMPACT OF THERMAL CORREC-TIONS
This section outlines the impact of T = 0 on the masses of particles in this model.
The formalism we follow is elaborately discussed in the review [33]. Henceforth, the thermal correction to the mass of a particle P will be denoted by δm 2 P (T ) and its thermally corrected mass by M P (T ). One then notes M P (T ) = m 2 P + δm 2 P (T ) where m P is the mass for T = 0. We first discuss the correction to the Z µτ mass. The contribution coming from a complex scalar carrying a charge q S (see left panel of Fig. 2) is given by We add here that only the longitudinal component of a gauge boson receives thermal corrections. Similarly, the contribution coming from a chiral fermion (see right panel of Fig. 2),  As detailed in the previous sections, the DM φ has feeble interactions with the scalars as well as with the gauge field Z µτ . The smallness of such interaction strengths implies that the interaction rate of DM remains smaller than the Hubble expansion rate throughout the thermal course of the universe. Consequently, the DM φ is never in equilibrium with the thermal bath and is injected into the thermal plasma via annihilations and decays of other particles. This is the freeze-in mechanism in a nutshell. Of these, the dominant contribution comes from the decay since the annihilations typically undergo suppressions by propagators and additional couplings. The impact of the annihilations is hence neglected in this study hereafter. Also, since the DM φ does not enter the thermal bath at any point in its cosmological history, it is cold. Therefore, thermal corrections to the DM mass itself become negligible in this setup.
We remind here that some previous studies [39] have looked at the freeze-in dynamics for the U (1) Lµ−Lτ model in detail. As mentioned earlier, we shall refer to the standard picture that has emerged from such studies as "standard freeze-in" or SFI. Our goal in this study is to demonstrate the deviation from SFI when thermal corrections to the masses of both the decaying particle and the DM are taken. We assume that all the decaying particles (see Fig. 3), i.e., the two scalars and Z µτ , are throughout in thermal equilibrium.
The Boltzmann equation predicting the DM yield is then given by Here M Pl denoting the temperature and the expansion rate of the universe respectively. In addition, Y φ = n φ s refers to the comoving number density of the DM with s being the entropy density. Y eq i signifies the equilibrium densities with i = Z µτ , S 1 , S 2 . Here, in theory, S 1 , S 2 generically denote the two neutral scalars. It is pointed out that {S 1 , S 2 } respectively coincide with {h, S} and {h 1 , h 2 } before and after the spontaneous breakdown of the gauge symmetry. The thermal decay width for A → B C reads where K n (x) is the nth order modified Bessel function. The late-time DM yield Y φ (x ∞ ) is calculated by solving the Boltzmann equation. The DM relic density is then obtained using We mention here we have neglected the effect of kinetic mixing (a T = 0 phenomenon) on freeze-in production as it is radiatively suppressed [44] and therefore small. It is reminded that the parameters controlling the interaction strengths of φ and ultimately Y φ (x ∞ ) are λ Hφ , λ Sφ , g µτ and n µτ . For a clearer understanding of the interplay of the different thermal that the lifetime corresponding to the Z µτ → ν µ ν µ , ν τ ν τ decays is smaller than the age of the Universe by several orders of magnitude for the said choice of g µτ . Thus Z µτ is not cosmologically stable and thus does not contribute to the relic density. Finally, we would also like to emphasize that the interaction strength of Z µτ with the SM fermions is large enough to keep it in equilibrium in the early universe.
A. λ Hφ , λ Sφ << n µτ g µτ This limit entails that Γ S i →φφ † << Γ Zµτ →φφ † and hence DM is dominantly produced by the Z µτ decay. In order to study the impact of thermal corrections on DM production, we propose m φ = 1 GeV, 30 MeV, and 25 MeV as benchmark values. The rationale behind choosing such values will become clear in the ensuing discussion.
The thermally corrected mass of Z µτ for T > v µτ reads M Zµτ (T ) 5 3 g µτ T since n µτ 1 for freeze-in. It is once again reminded that while the mass of Z µτ entirely comes from thermal corrections for T > v µτ , the DM φ does get a bare mass equalling µ φ from the scalar potential in the said temperature range. Fig. 4 shows the variation of M Zµτ (x) for the chosen values of m φ . At a very high temperature, say T initial = 10 5 GeV (x initial = 10 −5 for m φ = 1 GeV), one finds M Zµτ (T initial ) = 64.55 GeV. However, this mass gap diminishes with decreasing T (increasing x) and a crossover is observed at T = T cr 1 (x = x cr 1 ) obtainable through M Zµτ (T cr 1 ) = 2m φ , beyond which M Zµτ (T ) < 2m φ . We hereafter refer to this as the first crossover. Using the expressions for M Zµτ (T ) given above, one derives T cr This crossover is seen to happen for m φ = 1 GeV and 30 MeV. We next come to discuss the role of v µτ in this scenario. An inspection of Table III reveals that T cr 1 > v µτ = 80 GeV for m φ = 1 GeV and 30 MeV. The spontaneous breaking of U (1) Lµ−Lτ takes place at T = v µτ thereby generating a squared mass equalling g 2 µτ v 2 µτ for Z µτ . The thermally corrected mass for the same therefore shows a kink at T = v µτ , as can be seen in Fig. 4. And this kink opens up the possibility of having M Zµτ (T = v µτ ) > 2m φ for a second time during the thermal evolution of this scenario. We compute the Z µτ thermal mass at this symmetry-breaking threshold in Table III and discover that this indeed happens in case of m φ = 30 MeV. Moreover, at some T cr 2 < v µτ , one again might encounter M Zµτ (T cr 2 ) = 2m φ for a second time. This is referred to here as the second crossover with the corresponding temperature being T cr 2 = . We mention here again that all crossover possibilities and the corresponding temperatures and x-values are summarised in Table III for     First, we point out that in the absence of thermal corrections, Z µτ is either massless (for T > v µτ ) or at best has a 40 MeV mass (T < v µτ ). That is, m Zµτ ≤ 2m φ in either case for all the m φ values chosen. Therefore, the decay Z µτ → φφ † remains kinematically closed, and no DM production must take place in the SFI picture. However, as detailed before, this decay mode kinematically opens up upon incorporating the thermal corrections, and DM production gets triggered and continues up to x = x cr 1 for the DM masses permitting the first crossover. The decay threshold closes at x = x cr 1 , and DM production abruptly stops causing the DM yield to saturate at Y φ (x cr 1 ) immediately after. For m φ = 1 GeV, the Z µτ → φφ † threshold does not reopen at any later point. And this explains the horizontal line to the right of x cr 1 in Fig. 5. For m φ = 30 MeV, DM production stops at the first crossover point, thereby causing the plateau immediately to the right of x cr 1 = 3.22 × 10 −4 . For this BP, the Z µτ → φφ † threshold reopens shortly after at the symmetry breaking point, and freeze-in production kicks in again.
This reopening is what shows up as the kink in Fig. 6 around x = 3.75 × 10 −4 . However, this second phase of DM production is rather short-lived and terminates permanently at the second crossover point. Hence a second horizontal region x cr 2 = 4.33 × 10 −4 onwards in Fig.  6.
Lastly, we discuss the freeze-in production for m φ = 25 MeV. In this case, DM production is unimpeded up to the second crossover point. Therefore, one expectedly finds a plateau starting at x cr 2 = 5.37 × 10 −4 in Fig. 7. The symmetry breaking only leads to the minor cusp around the corresponding x-value, i.e., x = 3.12 × 10 −4 .
For the DM φ being generated through the decay of Z µτ only, Y φ (x ∞ ) and therefore Ω φ h 2 ∝ n 2 µτ . Having displayed the various temperature thresholds in the Figs. 5, 6 and 7, we subsequently tune n µτ appropriately such that Y φ (x ∞ ) is in the requisite ∼ 10 −11 ballpark.  In this section, we study the impact of the h 1 , h 2 → φφ † decays on freeze-in production for the chosen m φ values. It, therefore, becomes pertinent here to examine the thermal corrections to the masses of the decaying scalars. A crucial difference between Z µτ → φφ † and h 1 , h 2 → φφ † in this model is that while the former can lead to DM production at very early epochs (or at a very high T ) through thermal corrections, the latter is triggered primarily through spontaneous symmetry breaking. Given that we have v µτ = 80 GeV in addition to v = 246 GeV, one can treat v SB ∼ 100 GeV as a common symmetry breaking scale, and therefore the scalar decays are activated for T ≤ v SB . We further take M h 2 = 25 GeV, sinθ = 0.01 and y eµ = y eτ = 0.5 consistently with the collider constraints and neutrino data. And the impact of thermal corrections on the scalar masses becomes subdominant for such a temperature range. We first quote below the thermally corrected masses for h 1,2 to test this impact. Neglecting the effect of a small s θ , and the couplings λ Hφ and λ Sφ , one writes Such choices for m φ and m h 2 , therefore, imply that unlike Z µτ → φφ † , the h 1 , h 2 → φφ † decays can occur even in the absence of temperature effects. Therefore, DM production in the SFI picture is not completely ruled out for this subsection. greater than 2m φ in the T < v µτ range. In other words, the finite temperature corrections do not alter the original hierarchy in this case. We also take λ Hφ = λ Sφ ≡ λ for simplicity.
For m φ = 1 GeV, DM production from scalar decay kicks in around x 0.01, an epoch when the production from Z µτ decay has already ceased long back. As can be seen in Fig. 8, the scalar decays thus "lift" the horizontal DM yield curve, and the natural freeze-in saturation smoothly is attained around x ∼ O(0.1). No new crossovers are introduced in the process.
The cases of BP2 and BP3 are also not very different. For BP2 and BP3, DM production from Z µτ decay stops at x cr 2 = 4.33×10 −4 and 5.37×10 −4 respectively. However h 1 , h 2 → φφ † get activated at a slightly earlier epoch, i.e., x 3 × 10 −4 . So, DM matter production never entirely ceases for BP2 and BP3. This is corroborated by Figs. 9 and 10. Now that scalar decays too contribute to DM production, the parameters λ and n µτ values need to be chosen   carefully so as to obtain the required relic density. These parameters values for each BP can be read from the corresponding figure. For convenience, they are also listed in Table IV.
The same table also shows the percentage contribution to the observed relic through SFI for the three BPs. These numbers show that thermal correction to the mass of the decaying Z µτ can indeed lead to sizeable increments in the relic. Lastly, we state for completeness that the values listed for λ and n µτ in Table IV lead to direct detection cross sections well below the LZ limit.
Next, we use Fig. 11 to present a clear distinction between the SFI and FFI scenarios. We DM production are the SFI channels (resulting from the scalar decays). As expected, one needs to decrease the value of quartic coupling λ for an increasing DM mass in order to satisfy the observed DM relic density. The pink region above the red solid line shows the region of parameter space where the DM is always overproduced through the SFI channels. The purple and orange solid lines correspond to a parameter space where the correct DM relic density is observed but as a result of an interplay between SFI and FFI for n µτ = 1.3 × 10 −3 and n µτ = 1.46 × 10 −3 respectively. The light orange region shows to parameter space where the relic again remains overabundant. Finally, the different colored stars correspond to our BPs.
Whenever n µτ is non-zero, DM production from the decay of Z µτ through FFI also comes into the picture. Now, FFI together with the SFI helps the DM to satisfy the observed DM relic abundance. This behavior is also evident from Fig. 8, 9 and 10 which corresponds to our BP1, BP2 and BP3. These BPs are shown using different colored stars in Fig. 11. It is interesting to point out that increasing n µτ gradually increases the contribution of Z µτ decay towards the DM abundance and hence the contributions from the scalar decay (SFI) has to be appropriately tuned in order to satisfy the correct DM relic abundance. This behavior is depicted by the purple curve where we set n µτ = 1.3 × 10 −3 . We notice that this curve remains independent of λ for any λ ≤ 8 × 10 −12 for DM mass 24 MeV. This is because, for this combination of n µτ and m φ , the DM production from the Z µτ decay can itself give 100% contribution towards the DM relic density (production of DM from the scalar decay remains subdominant). Increasing λ above 8 × 10 −12 can results in an overproduction of the DM as the contribution from scalar decay (SFI channels) towards DM production also becomes significant. Once the peak is reached the curve starts to fall with the increasing DM mass. This is expected as in this regime, both FFI and SFI channels contribute towards the DM production to satisfy the correct relic abundance. Increasing n µτ further leads to an increase in the DM relic abundance for a fixed DM mass and hence a smaller λ is required in order to satisfy the observed DM relic abundance. This is clear from the orange curve where n µτ = 1.46 × 10 −3 . The trend observed for non-zero n µτ in m φ − λ plane also suggests that there exists an upper bound on the combination of n µτ and m φ (which corresponds to a vertical line) that can give the correct DM relic abundance with a dominant contribution coming from the decay of Z µτ and going beyond it always results in an overabundance of DM (light orange region). For this particular choice of fixed parameters, we found the combination to be n µτ = 1.46 × 10 −3 and m φ = 32 MeV (orange curve). Next, we would also like to point out the non-trivial shift observed in the DM mass when n µτ is increased from a smaller to a larger value (shift in the vertical line). Here, we notice that a larger n µτ requires a larger DM mass to satisfy the observed DM relic density. The reason for this behavior is the following, for a larger DM mass the first crossover of M Zµτ (T ) and 2m φ takes place at an earlier epoch in comparison to what happens for a smaller DM mass. This early crossover results in a smaller DM abundance when produced from the decay of Z µτ . Hence demanding the observed relic density accordingly requires a larger n µτ .
Finally, we will also like to comment briefly on the possibility of having a large v µτ with the same choice of g µτ fixed at 5 × 10 −4 as discussed above. A large v µτ results in a heavier Z µτ (vertical part of the orange curve). Although this scenario does not remain consistent with the recent g − 2 data, as can also be seen from Fig. 1, it still is interesting in terms of DM phenomenology. If the mass hierarchy among the Z µτ , the scalar responsible for breaking U (1) Lµ−Lτ symmetry, and the DM are appropriately set the scenario can result in the forbidden production of the DM from the decay of this scalar at a very early epoch through thermal corrections. The production of DM from Z µτ decay will always proceed through the SFI even with the thermal mass of Z µτ is taken into account if the DM mass is smaller than half of the Z µτ mass. Discussing this scenario in detail is beyond the scope of the present work, and we wish to take it as a future project.

V. SUMMARY AND CONCLUSION
Although the studies of FIMP dark matter in a minimally extended U (1) Lµ−Lτ model already exists in the literature, the role of thermal corrections in this setup has never been examined before. In this work, we show that incorporating a thermal mass for the gauge boson Z µτ opens up new temperature thresholds that are not encountered in SFI. For simplicity, we only consider the production of the DM through the decay of the gauge boson of U (1) Lµ−Lτ symmetry and two scalars. All the above-mentioned particles remain in equilibrium with the SM bath but couple feebly to the DM. The DM mass does not receive thermal corrections on account of the fact that a FIMP does not equilibrate with the thermal bath at any point. However, the masses of the decaying particles receive such corrections at high temperatures due to their interactions with the thermal bath.
For a better understanding of the role of different thermal masses, we divide our study into two different scenarios: (A) λ Hφ , λ Sφ << n µτ g µτ and (B) λ Hφ , λ Sφ ∼ n µτ g µτ . In the first scenario, the DM is dominantly produced by the decay of Z µτ whereas the second scenario entails the production of DM from the decay of Z µτ as well as the other two scalars. An exciting feature of the first scenario is the existence of two crossovers where the condition M Zµτ (T ) > 2m φ is satisfied. While the DM production always proceeds via a channel that remains kinematically forbidden in the SFI scenario before the first crossover, and the production of the DM after the second crossover might or might not happen via the forbidden channel. In the second scenario, the impact of the scalars decaying into the DM is also observed on top of its production from the decay of Z µτ . Here, the production from the scalar decay is triggered only after the spontaneous symmetry breaking. Finally, we also comment on the possibility of having a larger U (1) Lµ−Lτ breaking scale, which we wish to take as a future endeavor.
The involvement of the feeble interactions of DM with the SM particles in the freezein scenario makes the model exceedingly challenging to observe experimentally. With the WIMP DM parameter space almost getting ruled out from the present experiments, the FIMP-type DM has emerged as a new alternative. Hence, providing a detection prospect of such DM is always a compelling task. Keeping this in mind, we focused on a DM parameter of the present setup that remained consistent with the DM relic density but at the same time also provided a solution to the muon (g − 2) anomaly. This, in turn, also increases the predictability of the present setup.