Constraints on long-lived electrically charged massive particles from anomalous strong lens systems

TL;DR

This study uses strong lensing anomalies to constrain properties of long-lived electrically charged particles (CHAMPs), finding that their lifetime and abundance are limited by observed lensing effects, with implications for particle physics models.

Contribution

The paper develops a new method combining N-body simulations and lensing data to constrain CHAMP properties, a novel approach in this context.

Findings

For full CHAMP abundance, lifetime < 0.96 years at 95% CL.

Longer lifetimes are allowed if CHAMP fraction is less than 50%.

Constraints inform particle physics models involving CHAMPs.

Abstract

We investigate anomalous strong lens systems, particularly the effects of weak lensing by structures in the line of sight, in models with long-lived electrically charged massive particles (CHAMPs). In such models, matter density perturbations are suppressed through the acoustic damping and the flux ratio of lens systems are impacted, from which we can constrain the nature of CHAMPs. For this purpose, first we perform $N$-body simulations and develop a fitting formula to obtain non-linear matter power spectra in models where cold neutral dark matter and CHAMPs coexist in the early Universe. By using the observed anomalous quadruple lens samples, we obtained the constraints on the lifetime ($\tau_{\rm Ch}$) and the mass density fraction ($r_{\rm Ch}$) of CHAMPs. We show that, for $r_{\rm Ch}=1$, the lifetime is bounded as $\tau_{\rm Ch} < 0.96\,$yr (95% confidence level), while a longer lifetime $\tau_{\rm Ch} = 10\,$yr is allowed when $r_{\rm Ch} < 0.5$ at the 95% confidence level. Implications of our result for particle physics models are also discussed.

Tables

Table 1. Table 1 : Simulated models.

Model	${\tau }_{\mathrm{Ch}}$ [yr]	$r_{\mathrm{Ch}}$
CDM	-	0
CH(0.1, 0.1)	0.1	0.1
CH(1, 0.1)	1	0.1
CH(10, 0.1)	10	0.1
CH(0.1, 0.5)	0.1	0.5
CH(1, 0.5)	1	0.5
CH(10, 0.5)	10	0.5
CH(0.1, 0.7)	0.1	0.7
CH(1, 0.7)	1	0.7
CH(10, 0.7)	10	0.7
CH(0.1, 0.85)	0.1	0.85
CH(1, 0.85)	1	0.85
CH(10, 0.85)	10	0.85
CH(0.1, 1)	0.1	1
CH(1, 1)	1	1
CH(10, 1)	10	1

Table 2. Table 2 : Simulation setups. L 𝐿 L is the length on a side of the simulation box, N 𝑁 N is the number of the simulation particles, and ϵ italic-ϵ \epsilon is the gravitational softening length. HL10 is available only for the CDM, CH(1, 0.85), CH(10, 0.5), CH(0.1, 1), and CH(1, 1) models.

Setup	$L\left[\mathrm{Mpc}/h\right]$	$N$	$ϵ\left[\mathrm{kpc}/h\right]$
L5	5	${512}^3$	0.5
L10	10	${512}^3$	1.0
HL10	10	${1024}^3$	0.5

Equations

$$k_{\mathrm{Ch}}=aH{|}_{t={\tau }_{\mathrm{Ch}}}\simeq 0.18/Mpc{\left(\frac{1\mathrm{yr}}{{\tau }_{\mathrm{Ch}}}\right)}^{1/2},$$

$$k_{\rm J}=a\sqrt{\frac{4\pi G\rho_{\rm M}}{\sigma^{2}}}\bigg{|}_{t=t_{\rm eq}}=20/{\rm Mpc}\,\left(\frac{m_{\rm WDM}}{1\,{\rm keV}}\right)^{4/3}\,,$$

$${\tau }_{\mathrm{Ch}}\simeq 1.5\mathrm{yr}{\left(\frac{1\mathrm{keV}}{m_{\mathrm{WDM}}}\right)}^{8/3}.$$

$$T^2(k_{\mathrm{d}})=\frac{1}{(1+k/k_{\mathrm{d}}{)}^{0.7441}},k_{\mathrm{d}}({\tau }_{\mathrm{Ch}})=503.2h/\mathrm{Mpc}{\left(\frac{1\mathrm{yr}}{{\tau }_{\mathrm{Ch}}}\right)}^{0.9273}.$$

$$\sum _{\begin{array}{c}{r}_{\mathrm{Ch}}=1\\ \mathrm{models}\end{array}}\sum _{z\mathrm{bins}}\sum _{k\mathrm{bins}}{\left(T^2-P_{\mathrm{Ch}}/P_{\mathrm{CDM}}{|}_{L5}\right)}^2,$$

$$T^2(f_{\mathrm{Ch}},k_{\mathrm{d}})=(1-f_{\mathrm{Ch}})+\frac{f_{\mathrm{Ch}}}{(1+k/k_{\mathrm{d}}{)}^{0.7441}},f_{\mathrm{Ch}}(r_{\mathrm{Ch}})=1-\exp \left(-1.382\frac{r_{\mathrm{Ch}}^{0.5168}}{1-r_{\mathrm{Ch}}^{0.9417}}\right),$$

$$\eta\equiv\biggl{[}\frac{1}{2N_{\rm pair}}\sum_{i\neq j}\left[\delta^{\mu}_{i}(minimum)-\delta^{\mu}_{j}(saddle)\right]^{2}\biggr{]}^{1/2}\,,$$

$$p(\tau_{\rm Ch},r_{\rm Ch})=\left(\int_{{\hat{V}}}\prod_{a}d\eta_{a}P(\eta_{a};\langle\eta_{a}^{2}\rangle^{1/2},\delta{\hat{\eta}}_{a})\right){\bigg{/}}\left(\int\prod_{a}d\eta_{a}P(\eta_{a};\langle\eta_{a}^{2}\rangle^{1/2},\delta{\hat{\eta}}_{a})\right)\,,$$

$$\prod _{a}P({\eta }_a\in \hat{V};⟨{\eta }_a^2{⟩}^{1/2},\delta {\hat{\eta }}_a)<\prod _{a}P({\hat{\eta }}_a;⟨{\eta }_a^2{⟩}^{1/2},\delta {\hat{\eta }}_a),$$

$$P(\eta ;⟨{\eta }^2{⟩}^{1/2})\propto \exp \left[-\left\{\ln (1+\eta /{\eta }_0(⟨{\eta }^2{⟩}^{1/2}))-\ln (\mu )\right\}/(2{\sigma }^2)\right]/(\eta +{\eta }_0).$$

$${\eta }_0(⟨{\eta }^2{⟩}^{1/2})=0.228⟨{\eta }^2{⟩}^{1/2},\mu =4.10,{\sigma }^2=0.279,$$

$$\langle\eta^{2}\rangle=\frac{1}{4}\biggl{[}(I_{A}+I_{B})-2I_{AB}(\theta_{AB})+(I_{B}+I_{C})-2I_{BC}(\theta_{BC})\biggr{]}\,,$$

$$I_{ij}({\theta }_{ij})$$

$$\langle\delta\kappa(0)\delta\kappa(\theta)\rangle=\frac{9H_{0}^{4}\Omega_{m,0}^{2}}{4c^{4}}\int_{0}^{r_{S}}drr^{2}\biggl{(}\frac{r-r_{\rm S}}{r_{\rm S}}\biggr{)}^{2}[1+z(r)]^{2}\int_{k_{\rm{lens}}}^{k_{\textrm{max}}}\frac{dk}{2\pi}kW_{\textrm{CS}}^{2}(k;k_{\textrm{cut}})P_{\delta}(k,r)J_{0}(g(r)k\theta)\,,$$

$$g(r)=\left\{\begin{array}[]{ll}r&\mbox{for $r<r_{\rm L}$}\\ {r_{L}(r_{S}-r)}/{(r_{S}-r_{L})}&\mbox{for $r\geq r_{\rm L}$}\end{array}\,,\right.$$

Full text

1 Introduction

Long-lived electrically charged massive particles (CHAMPs) [1, 2, 3, 4] are suggested to resolve the failure of the standard $\Lambda$CDM ($\Lambda$: dark energy; CDM: cold dark matter) model in reproducing the observed number of satellite galaxies in the Milky Way [5, 6] with satisfying cosmological bounds, e.g., from big-bang nucleosynthesis [7, 8, 9, 10, 11, 12], cosmic microwave background (CMB) anisotropies [13], and magnetic fields in galaxy clusters [14] (for other bounds, see [13] and references therein). The Coulomb interaction ties CHAMPs to baryons until CHAMPs decay into neutral, stable, and massive particles, which account for part of the present DM mass density. CHAMPs are prevented from falling into the bottom of the primordial gravitational potential by the induced pressure. The resultant matter power spectra are suppressed below the comoving horizon length scale at the CHAMP decay when compared to that in the standard CDM model.

As a concrete example of long-lived CHAMPs in particle physics models, especially when we consider supersymmetric (SUSY) models, it is known that a scalar partner of charged leptons tends to have a electroweak-scale mass and a long lifetime. Two possible cases have been known, where the scalar charged leptons can indeed become long-lived in the minimal-supersymmetric standard model (MSSM) (see, e.g., [15, 16, 17, 18]). First, the lightest SUSY particle (LSP) can be the gravitino with the next-to-LSP (NLSP) being a scalar charged lepton, e.g., the scalar tau lepton (i.e., stau) in the gauge-mediated SUSY breaking models [19, 20]. The NLSP stau decays into the gravitino LSP only through the gravitational interaction and thus the electroweak-scale stau has a long lifetime. Second, in the case of a degenerate mass spectrum between the stau ($\tilde{\tau}$) NLSP and the neutralino ($\chi$) LSP, the stau can be long-lived kinematically due to the smallness of the mass difference. In the latter model, it is remarkable that the lifetime can be of the order of $0.1\text{--}10$ yr with their small mass difference being $\Delta m<50$ MeV and their large absolute masses being $\gtrsim 100$ GeV [21]. In this case, a possible decay mode should be only into $\tilde{\tau}\to\chi+e+\nu_{e}+\nu_{\tau}$ with a negligible kinetic energy of the daughter electron ($e$) [21]. Therefore, the decay of the stau NLSP does not induce any significant photodissociations of background 4He or D [22, 23] at a low temperature of $T\lesssim 10$ keV, which agrees with the observed light element abundances.

Even if we consider a more general case that is no longer based on the MSSM, a degenerate mass spectrum between the CHAMP and its daughter particle kinematically induces a long lifetime of the CHAMP, while neither high-energy charged leptons nor photons are emitted through the decay by a small mass difference. Remarkably, the lifetime of the CHAMP tends to be longer for a larger mass of the CHAMP with the mass difference being kept fixed. For example, in order to obtain a conservative bound on the CHAMP, it would be possible to consider a mass spectrum in which $\Delta m$ is smaller than ${\cal O}(10)$ MeV by which any photodissociations of the light elements are not induced, and the free-streaming length of the daughter particle is too short to erase the density perturbations below $k={\mathcal{O}}(10^{4})\,h/$Mpc. In this paper therefore, we consider such conservative setups for the decaying CHAMP models by assuming that negligible amounts of energetic electromagnetic particles are emitted with a negligible kinetic energy of the daughter particle.

We use the analyzing method in [24, 25] though we consider long-lived CHAMP models instead of warm (WDM) or mixed (MDM) dark matter models. In both the models the formation of subgalactic-scale ($k={\mathcal{O}}(100)\,h$/Mpc) objects is suppressed. One may wonder if we can map a derived constraint on MDM models, more specifically, the WDM mass and the mass density fraction of the warm component ($m_{\rm WDM}$ and $r_{\rm warm}$) to the lifetime and the mass density fraction of CHAMPs ($\tau_{\rm Ch}$ and $r_{\rm Ch}$) with a help of the relation derived in WDM models [4]: $\tau_{\rm Ch}\simeq 1.5\,{\rm yr}\,(1\,{\rm keV}/m_{\rm WDM})^{8/3}$. Our result actually confirms that this relation is valid between pure long-lived CHAMP models ($r_{\rm Ch}=1$) and WDM models ($r_{\rm warm}=1$). Nevertheless, as we will see later, this mapping is not applicable between mixed long-lived CHAMP ($r_{\rm Ch}<1$) models and MDM ($r_{\rm warm}<1$) models. This is virtually because of the difference in suppression mechanisms between long-lived CHAMP models and MDM models: the acoustic damping through the (indirect) interaction with relativistic species and the collisionless damping through the free-streaming of the warm component, respectively. Thus we perform independent analysis in long-lived CHAMP models in this paper, even though the methodology itself is close to that in [24, 25].

In order to constrain long-lived CHAMP models, we use quadruply lensed quasar-galaxy systems that show anomalous flux ratios in lensed images; the flux ratios of lensed images disagree with the prediction of the best-fit lens models with a smooth potential whose variation length scale is larger than the separation between the lensed images. It has been argued that such anomalies are caused by subhalos in lens galaxies [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. However, recent analyses showed that the contribution from intervening halos are significant [39, 40, 41, 42], while the predicted subhalo population is too low to explain the observed anomalous flux ratios [43, 44, 36, 37, 45, 46]. The contribution from subhalos is turned out to be just $\sim 30\%$ of the total that includes contribution from intervening structures such as halos, voids, and filaments along the line of sight [47, 48]. In fact, a sign of perturbations by locally negative density perturbations has been observed in the quadruply lensed submillimeter galaxy (SDP.81) [49]. In our analysis, we assume that the observed anomalous flux ratios are mainly caused by intervening structures rather than subhalos in lens galaxies [50, 41].

We use optical/near-infrared data for the positions of lensed images and the centroid of the foreground elliptical galaxy and mid-infrared/radio (cm) data for the flux of lensed images. Since the angular size of lensed images in the mid-infrared/radio bands is sufficiently larger than the Einstein angular radius of stars in the foreground galaxy, the change in the flux of lensed images due to microlenses is negligible in the mid-infrared/radio bands. We assume that the primary lens in each system is described by a smooth potential with an elliptical symmetry in the surface mass density projected onto the lens plane. We fit the data with a model that includes a primary lens with an external shear and possibly a companion galaxy if necessary. Then the magnification perturbation in each system is estimated from the residual in the flux of lensed images obtained from the fit. If the acoustic damping of small-scale density perturbations is too strong, then the magnification perturbation cannot be explained by structures in the line of sight. Thus we can constrain the lifetime and the mass density fraction of CHAMPs.

The organization of this paper is as follows. In section 2 we study the behavior of density perturbations in both the linear and non-linear regimes. There we describe the physical mechanism that results in the suppression of the linear matter power spectrum in long-lived CHAMP models. A fitting formula of non-linear matter power spectra is derived. By using the fitting formula we evaluate the second moment of the probability density function (PDF) of the magnification perturbation in section 3. Then we derive the probability of reproducing the observed flux ratios in a given long-lived CHAMP model, which is used to constrain the long-lived CHAMP model parameters. Section 4 is devoted to concluding remarks. Throughout this paper, we take the cosmological parameters obtained from the observed CMB anisotropies (“Planck + WP” (WP: WMAP polarization) [51]) to be consistent with [24, 25]: the matter mass density at present, $\Omega_{m,0}=0.3134$; the baryon mass density, $\Omega_{b,0}=0.0487$; the cosmological constant, $\Omega_{\Lambda,0}=0.6866$; the Hubble constant, $H_{0}=67.3$ $(=100h)\,{\rm km/s/Mpc}$; the spectral index, $n_{s}=0.9603$; the root-mean-square (rms) amplitude of matter density perturbations at $8h^{-1}\,{\rm Mpc}$, $\sigma_{8}=0.8421$.

2 Matter power spectra

2.1 Linear matter power spectra

We suitably modify the public code CAMB [52] to incorporate the evolution equations of the long-lived CHAMP perturbations given in [1, 2]. Here the tight coupling between CHAMPs and baryons is assumed, which is valid as long as $m_{\rm Ch}<10^{8}$ GeV [13]. In this mass range, the matter power spectra in long-lived CHAMP models are parametrized by the lifetime and the mass density fraction: $\tau_{\rm Ch}$ and $r_{\rm Ch}=\Omega_{{\rm Ch},0}/\Omega_{c,0}$, where $\Omega_{{\rm Ch},0}$ is the CHAMP mass density evaluated assuming that it is stable until today even for the case with the lifetime ($\tau_{\rm Ch}$) is shorter than the age of the Universe. $\Omega_{c,0}$ is the CDM mass density at present. We assume that some CDM-like particles account for the rest of the observed DM mass density: $1-r_{\rm Ch}$. Note that this CDM-like particle could be identical to the neutral decay product of the CHAMP, but their origins and thermal histories are completely different; the CDM-like particle should exist as a stable and cold component of the Universe at least around and after the horizon entry of the perturbations of interest. A cutoff scale can be estimated at the comoving horizon scale at the CHAMP decay:

[TABLE]

with the scale factor $a$ being normalized such that $a=1$ at present and the Hubble expansion rate being $H$.

In general $r_{\rm Ch}$ can be larger than unity with the present mass density of the neutral decay product equal to or smaller than the observed DM mass density: $\Omega_{{\rm daughter},0}/\Omega_{c,0}\leq 1$. This is because the CHAMP mass can be much larger than the mass of the neutral CHAMP decay product. In this paper, on the other hand, we focus on the case that the two masses are degenerate: $m_{\rm Ch}\simeq m_{\rm daughter}$. This choice has two advantages: natural explanation of the CHAMP long lifetime and harmlessness of charged decay product. One may wonder how we can realize a relatively long-lived massive particle. One possible reason is that the electroweak-scale CHAMP decays into the neutral decay product through a superweak interaction, say, gravitationally (see [15, 16, 17, 18] for supersymmetric realizations). When the CHAMP mass is much larger than the electroweak scale, even a gravitational decay is not sufficient to make the CHAMP stable over ${\mathcal{O}}(1)$ yr. The small mass difference between the CHAMP and the neutral decay product restricts the final-state phase space and lengthens the CHAMP lifetime. We also care about the energy injection from the CHAMP decay into the standard model plasma, which may be disfavored since they could destroy the light element produced from the big-bang nucleosynthesis [22, 23, 12], and/or distort the CMB black body spectrum to generate the $y$- and $\mu$-distortions [53, 54].

Figure 1 shows the linear matter power spectra extrapolated to the present time ($z=0$) in some long-lived CHAMP, WDM, and MDM models. For notational simplicity, we denote a long-lived CHAMP model with the lifetime of $\tau_{\rm Ch}$ and the mass density fraction of $r_{\rm Ch}$ as CH($\tau_{\rm Ch}$,  $r_{\rm Ch}$), where $\tau_{\rm Ch}$ is measured in units of year as shown in Table 1. For reference we also show the linear matter power spectra in MDM models, where the warm component follows the Fermi–Dirac distribution with the spin being one half and is characterized by the warm particle mass and the mass density fraction: $m_{\rm WDM}$ and $r_{\rm warm}=\Omega_{{\rm warm},0}/\Omega_{c,0}$. A cutoff scale in WDM models is given by the Jeans scale at matter-radiation equality:

[TABLE]

with the velocity dispersion being $\sigma^{2}$ [4]. The face values of eqs. (2.1) and (2.2) do not completely coincide with the cutoff scales read from calculated linear matter power spectra. Nevertheless the scaling with the parameter is found appropriate. It is suggested that we use $k_{\rm cut}\simeq 12k_{\rm Ch}\simeq k_{\rm J}/12$ as an empirical cutoff scale [4]. The mapping of the model parameters between pure long-lived CHAMP models and WDM models can be done with the following relation:

[TABLE]

As seen from figure 1, the CH(1, 1) model (orange) shows an oscillation, which is not seen in the WDM model with $m_{\rm WDM}=1.3$  keV (blue); the density perturbation crosses zero around $k\simeq 14\,h$/Mpc and the peak at $k\simeq 20\,h$/Mpc has the opposite sign to the perturbations at $k<10\,h$/Mpc. The oscillation is an imprint of the acoustic oscillation between CHAMPs and baryons that takes place until the CHAMP decay.

Let us highlight another important feature by comparing the CH(10, 0.5) model (green) and the MDM model with $(m_{\rm WDM}\,[{\rm keV}],r_{\rm warm})=(0.5,0.5)$ (red). The suppression at smaller length scales is weaker in the long-lived CHAMP models than in the MDM models even with the mass density fraction being identical: $r_{\rm Ch}=r_{\rm warm}$. This is because the suppression mechanisms are different between them. Long-lived CHAMPs suppress the density perturbations of the neutral decay product through the acoustic damping discussed above, while in MDM models a collisionless damping is essential and the free-streaming of the warm component transfers the power of the density perturbation into the higher order terms of the Boltzmann hierarchy. The former works in the radiation dominated era, or precisely speaking, until baryons including CHAMPs dominate the entropy density or the long-lived CHAMP decays, while the latter works both in the radiation dominated era and in the matter dominated era. In the radiation dominated era the gravitational coupling between long-lived CHAMPs and the cold component is negligible since the gravitational potential is determined by the radiation energy density. On the other hand, in MDM models, the gravitational force between the warm and cold components becomes important in the matter dominated era; the velocity dispersion of the warm component prevents not only the warm component but also the cold one from growing gravitationally. This observation is supported by figure 1, which shows that the matter powers in the CH(10, 0.5) model is a quarter $(=r_{\rm Ch}^{2})$ smaller than the CDM (black) at smaller length scales.

2.2 Non-linear matter power spectra

We perform $N$-body simulations by using linear matter power spectra calculated in the models listed in table 1 to generate the initial conditions. The public code Gadget-2 [55] is used with the initial redshift of $z=49$ and the other setups summarized in table 2. The simulation particles represent the whole matter distribution and are random samplings of its initial phase space distribution through the Zel’dovich approximation. The realization of the random variable is identical among the simulated long-lived CHAMP models as well as the simulated MDM models investigated in [25].

The direct use of the simulated matter distributions is preferable to calculate the magnification perturbations in long-lived CHAMP models. The sampling given in table 1, on the other hand, is not sufficiently dense to obtain a smooth constraint map in the $\tau_{\rm Ch}$-$r_{\rm Ch}$ plane. Relying on the weak lens approximation, we can evaluate the magnification perturbations from the non-linear matter power spectra. To this end we construct the fitting formula of the non-linear matter power spectra as a function of $\tau_{\rm Ch}$ and $r_{\rm Ch}$ through the following steps: we construct the fitting formula of the L5 and L10 simulated transfer functions ($T^{2}=P_{\rm Ch}/P_{\rm CDM}$); we obtain the fitting formula of the non-linear matter power spectra by multiplying the halofit [56] and $T^{2}$; and then we check the derived fitting formula with the HL10 simulations.

First we measure the power spectra of the simulated whole matter distributions in the L5 and L10 simulations in the pure long-lived CHAMP models ($r_{\rm Ch}=1$) at $z=0,0.3,0.6,1,2$, and $3$. Figure 2 shows the simulated transfer functions in the CH(0.1, 1) and CH(1, 1) models: $P_{\rm Ch}/P_{\rm CDM}|_{\rm L5,L10}$. The L10 simulations show slightly larger matter powers than the L5 simulations. This is because a smaller box size simulation generically misses effects of non-linear coupling between smaller and larger length-scale perturbations, which are more important in suppressed linear matter power models than in the standard CDM model. On the other hand, the agreement between the L5 and L10 simulations is up to 10% in the amplitude of the transfer function above the Nyquist wavenumber of the L5 simulation ($k_{\rm Nyq}=322\,h/$Mpc) at $z<2$ and sufficient for our purpose.

Interestingly the following fitting formula motivated by analogy to that in WDM models [24] can reproduce the pure long-lived CHAMP simulations up to 10% in the amplitude of the transfer function as shown in figure 2:

[TABLE]

Here we vary and determine the overall factor and the scaling index of $k_{\rm d}$, which is a function of $\tau_{\rm Ch}$, by minimizing the residual:

[TABLE]

where $z\,{\rm bins}$, $z\in\{0,0.3,0.6,1,2,3\}$, and $k\,{\rm bins}$, $\log(k\,[h/{\rm Mpc}])\in\{2.117,2.137,2.157,\cdots,2.477\}$ (total 20 bins), have a large number of samples and are not affected by cosmic variance. Note that we can evaluate the overall factor and scaling index as $543.8$ and $0.8277$, respectively, by using the relation given in eq. (2.3) though it is read from the linear matter power spectra. These numbers are quite close to those obtained by the above fitting procedure, which support the suggestion of [4]: the relation given in eq. (2.3) map WDM models to pure long-lived CHAMP models even after the non-linear growth.

Next we extend the above fitting formula to the transfer functions in mixed long-lived CHAMP models. We have found that the following extension by analogy to that from WDM models to MDM models [25] provides a good fit at the level of $\sim 10$% in amplitude:

[TABLE]

where $k_{\rm d}$ is the same as given in eq. (2.4). We assume a functional form of $f_{\rm Ch}$ as $f_{\rm Ch}(r_{\rm Ch})=1-\exp[-ar_{\rm Ch}^{b}/(1-r_{\rm Ch}^{c})]$ by analogy to $f_{\rm warm}$ in [25], and then vary and determine the three parameters ($a,b,$ and $c$) by minimizing the residual given in eq. (2.5) with the sum over the mixed long-lived CHAMP models ($r_{\rm Ch}<1$) listed in table 1.

The comparison between $T^{2}$ and $P_{\rm MDM}/P_{\rm CDM}|_{\rm L5,L10}$ is given in figure 3. Note that the MDM model with $(m_{\rm WDM}\,[{\rm keV}],r_{\rm warm})=(1,0.8)$ (blue), which shows the comparable transfer function to that in the CH(1, 0.85) model, is disfavored at the 95% confidence level. Thus we expect that the CH(1, 0.85) model is also excluded at the 95% confidence level, which provides a consistency check for our analysis below.

We construct the fitting formula of the non-linear matter power spectra in mixed long-lived CHAMP models by multiplying $T^{2}$ (eq. (2.6)) and halofit, which has been developed for the standard CDM model with the same cosmological parameters as employed in this paper. The explicit expression of halofit can be found in [58, 24] and thus is not repeated here. The agreement between halofit and simulations with our setups is checked and confirmed in [25]. We compare the fitting formula of the non-linear matter power spectra and the simulations in figures 4 and 5. Here we further check the fitting formula by using the LH10 simulations. The agreement between HL10 and L10 sumulations is within at worst 10% level in the non-linear matter powers, while L5 is 10% smaller than these two. A net error of the fitting formula can be estimated at 20% at most in amplitude at $z<2$.

3 Lens analysis and result

The anomalous quadruple lens samples and the analyzing method are the same as in [24, 25]. The samples are B1422$+$231, B0128$+$437, MG0414$+$0534, and B0712$+$472. The primary lens galaxy halo is modeled by a singular isothermal ellipsoid. The parameters in the model are determined by fitting the positions of lensed images and the centroid of the primary lensing galaxy. The derived flux ratio does not fit the observed data at the 95% confidence level or more.

To quantify the probability of reproducing the observed flux ratios in a given model, we use the method derived in [41, 42]. Below we briefly describe the methodology of our analysis. For details, we refer the readers to [41, 42]. The magnification perturbation can be characterized by an estimator of

[TABLE]

in which the sum is taken over $N_{\rm pair}$ pairs of images of a point source with different parities and $\delta^{\mu}_{i}=\delta\mu_{i}/\mu_{i}$ with $\mu_{i}$ and $\delta\mu_{i}$ being the magnification and its perturbation of the image $i$, which is evaluated at minimum and saddle points in the arrival time surface.

To obtain constraints on the CHAMP model parameters ($\tau_{\rm Ch}$ and $r_{\rm Ch}$), we calculate the $p$-value with the following formula:

[TABLE]

where $P(\eta_{a};\langle\eta_{a}^{2}\rangle^{1/2},\delta{\hat{\eta}}_{a})$ is the PDF for the anomalous lens system $a$ with $\eta_{a}$ and $\langle\eta_{a}^{2}\rangle$ being respectively the estimator for the system $a$ and the ensemble average of its second moment. $\delta{\hat{\eta}}$ is the observational error which is taken into account by replacing $\langle\eta^{2}\rangle^{1/2}$ as $\langle\eta^{2}\rangle^{1/2}\to(\langle\eta^{2}\rangle+\delta{\hat{\eta}}^{2})^{1/2}$ in the PDF. We define the integration domain (${\hat{V}}$) by

[TABLE]

with the observed value of the estimator being ${\hat{\eta}}$. For the PDF of $\eta$, we assume a lognormal form as follows:

[TABLE]

Here three parameters $\eta_{0},\mu$, and $\sigma$ are taken to be

[TABLE]

which have been calibrated to the PDF obtained by ray-tracing simulations in the standard CDM model [42]. We adopt the values given in eq. (3.5), which is considered to be a reasonable assumption in WDM models [24], also in long-lived CHAMP models.

Here we adopt the weak lens approximation, where $|\delta_{i}^{\mu}|\ll 1$. For the case with three images A, B, and C where A and B are minimum points and C is a saddle one, the second moment of the estimator can be calculated as

[TABLE]

where $\theta_{ij}$ $(i,j=\textrm{A,B,C})$ corresponds to the separation angles between the images $i$ and $j$. Here $I_{ij}(\theta_{ij})$ is defined by

[TABLE]

where $\kappa_{i}$ is the convergence, $\gamma_{1i}$ and $\gamma_{2i}$ are the components of the shear, and $\delta\kappa_{i},\delta\gamma_{1i}$, and $\delta\gamma_{2i}$ are their perturbations of the image $i$. Here we take the coordinates such that $\langle\delta\kappa\delta\gamma_{2}\rangle$ and $\langle\delta\gamma_{1}\delta\gamma_{2}\rangle$ vanish by setting the separation angle between a pair of lensed images perpendicular to the $+$ mode. $I_{i}$ $(i=\textrm{A,B,C})$ is defined by using eq. (3.7) as $I_{i}=I_{ii}(0)$.

The non-linear matter power spectrum ($P_{\delta}(k,r)$) at the comoving distance of $r$ is related to the auto-correlation function of $\delta\kappa$ as follows:

[TABLE]

where $z(r)$ is the redshift to the comoving distance and $J_{n}$ is the $n$th order Bessel function. $g(r)$ is given by

[TABLE]

which describes the photon trajectory in a primary lens, where $r_{L}$ and $r_{S}$ represent respectively the comoving distances to the lens galaxy and the source. In the $k$-integral of eq. (3.8), the lower bound is given by $k_{\textrm{lens}}=\pi/(2r_{\rm L}b)$ with $b$ being the mean angular separation between a lensed image and a lens centre. The upper bound is chosen to be $k_{\textrm{max}}=\mathcal{O}(10^{3\text{--}5})\,h/$Mpc, which corresponds to the scale above which perturbations become negligible due to the finite source size ($\sim\mathcal{O}(1)$ pc). $W_{\rm CS}$ is a filtering function which is called the constant shift (CS) filter [42] and mildly cuts the contribution from large angular-scale modes with $k<k_{\rm cut}$. Other correlation functions $\langle\delta\gamma_{1}(0)\delta\gamma_{1}(\theta)\rangle$, $\langle\delta\gamma_{2}(0)\delta\gamma_{2}(\theta)\rangle$, and $\langle\delta\kappa(0)\delta\gamma_{1}(\theta)\rangle$ are also given by similar expressions whose explicit forms can be found in [24] and thus are not repeated here.

We compute the non-linear matter power spectra (halofit multiplied by $T^{2}$ (eq. (2.6))) in 25 models including the long-lived CHAMP (total 15) models listed in table 1. We calculate $\langle\eta_{a}^{2}\rangle^{1/2}(\tau_{\rm Ch},r_{\rm Ch})$ for each lens system $a$ from eq. (3.6). We interpolate the resultant $p$-values linearly in the $(\tau_{\rm Ch},r_{\rm Ch})$ plane and show the result in figure 6. The pure long-lived CHAMP models with $\tau_{\rm Ch}>0.96$ yr fail in reproducing the observed flux ratios, $p(\tau_{\rm Ch},r_{\rm Ch}=1)<0.05$. The constraint on $\tau_{\rm Ch}$ becomes weaker for smaller $r_{\rm Ch}$. This is because for a given $\tau_{\rm Ch}$, the suppression of matter power spectra is milder for smaller $r_{\rm Ch}$. The CH(1, 0.85) model is excluded as expected from its similar non-linear matter power spectrum to that in the MDM model with $(m_{\rm WDM}\,[{\rm keV}],r_{\rm warm})=(1,0.8)$ (see discussion in section 2.2). The mixed long-lived CHAMP models with the lifetime as long as $\tau_{\rm Ch}=10$ yr are compatible if $r_{\rm Ch}<0.5$ at the 95% confidence level.

4 Conclusion

We have examined the possibility of reproducing the observed flux ratios in four anomalous quadruple lens systems via weak lens effects of the line-of-sight matter distribution in long-lived CHAMP models. In long-lived CHAMP models the density perturbations at subgalactic scales are suppressed as in WDM or MDM models, whose behavior of the magnification perturbations was studied in [24, 25]. Although our method for constraining CHAMP models is similar to that used for constraining MDM models, an independent analysis was required. In the two models, the suppression mechanisms of the matter power spectra are different: the acoustic damping in long-lived CHAMP models and the collisionless damping in MDM models. In fact, we can map WDM models to pure long-lived CHAMP models as suggested in the previous literature, but the mapping of the model parameters between mixed long-lived CHAMP models and MDM models is non-trivial. This point has been clearly seen in our calculated linear matter power spectra. Meanwhile the mapping of the model parameters between pure long-lived CHAMP models and WDM models has been found valid in the non-linear matter power spectra, or in other words, even after the non-linear gravitational evolution.

Our analysis has shown that the pure long-lived CHAMP models with $\tau_{\rm Ch}<0.96$ yr are successful in reproducing the observed flux ratios at the 95% confidence level or more. A shorter CHAMP lifetime implies an earlier realization of the CDM model and thus a wider range matches of the matter power spectra between long-lived CHAMP models and the standard CDM model. An upper bound on the CHAMP lifetime corresponds to a lower bound on the mass of the warm component in WDM models, where the CDM model is the heavy mass limit of the WDM model. The mapping of the model parameters between them is compatible with the previous literature: $\tau_{\rm Ch}\simeq 1.5\,{\rm yr}\,(1\,{\rm keV}/m_{\rm WDM})^{8/3}$.

Once the mass density fraction of the long-lived CHAMP component is reduced and some CDM-like particle is introduced to account for the rest, the constraint on the CHAMP lifetime has become weaker. This is because the suppressions in linear matter power spectra become $r_{\rm Ch}^{2}$ times milder when compared to pure long-lived CHAMP models ($r_{\rm Ch}=1$) with the same CHAMP lifetime. The derived constraints on the long-lived CHAMP mass density fraction is $r_{\rm Ch}<0.5$ (95% confidence level) for the CHAMP lifetime of $\tau_{\rm Ch}=10$ yr.

In order to obtain such a long lifetime to be sensitive to the current analyses in concrete models, e.g., in the MSSM, there are two possible scenarios of the stau NLSP decaying into a LSP dark matter particle: 1) into the gravitino LSP and 2) into the neutralino LSP with a small mass difference between the stau and the neutralino. In either scenario, the masses should be commonly degenerate between the CHAMP and the daughter particle to be $m_{\rm Ch}\sim m_{\rm daughter}$, which is a general condition to avoid the photodissociation of light elements [22, 23, 12] and the CMB distortion [53, 54], independently of particle physics models. Regarding the CHAMP abundance, when we adopt a non-thermal production of CHAMPs, we do not have to stick to the standard thermal relic abundance of CHAMPs. However, we have to be careful for the condition that the mass density of daughter particles produced by decays of CHAMPs should not exceed the observed DM mass density, i.e., $\Omega_{\rm daughter}\sim\Omega_{\rm Ch}=r_{\rm Ch}\Omega_{c,0}$. Meanwhile, to evade the constraint from the catalyzed effect on an overproduction of the 6Li abundance [7, 8, 9, 10, 11, 12], we have a strict big-bang nucleosynthesis bound on the ratio of the number density of CHAMPs to the entropy density ($s$): $n_{\rm Ch}/s\lesssim 10^{-15}$ for $\tau_{\rm Ch}\gtrsim 10^{3}$ s. In summary, we have a lower bound on the mass of the CHAMP, $m_{\rm Ch}\gtrsim 4\times 10^{2}r_{\rm Ch}$ TeV, which can be realized in concrete models with $r_{\rm Ch}<1$, e.g., in the MSSM.

Although, we have used only four samples of quadruple lens systems with flux-ratio anomalies, the number of such systems will surge in near future. For instance, the recently found quadruple lens systems of submillimeter galaxies [59] and Lyman-alpha emitting galaxies [60] are expected to contain many systems with flux-ratio anomalies. Using such new samples, we would have a more stringent constraint on the long-lived CHAMP model parameters.

Acknowledgments

This work was partially supported by IBS under the project code IBS-R018-D1 (AK), JSPS KAKENHI Grant Numbers 26247042 (KK) and 15K05084 (TT), and MEXT KAKENHI Grant Numbers 15H05889, 16H00877 (KK), and 15H05888 (TT). Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan and on SR16000 at YITP in Kyoto University.