Displaced new physics at colliders and the early universe before its first second

Displaced vertices at colliders, arising from the production and decay of long-lived particles, probe dark matter candidates produced via freeze-in. If one assumes a standard cosmological history, these decays happen inside the detector only if the dark matter is very light because of the relic density constraint. Here, we argue how displaced events could very well point to freeze-in within a non-standard early universe history. Focusing on the cosmology of inflationary reheating, we explore the interplay between the reheating temperature and collider signatures for minimal freeze-in scenarios. Observing displaced events at the LHC would allow to set an upper bound on the reheating temperature and, in general, to gather indirect information on the early history of the universe.


Introduction
Detecting non-gravitational interactions of dark matter (DM) has been the target of a vast and diverse experimental effort with no conclusive evidence so far [1][2][3]. In this work, we exploit the interplay between two different DM probes: the measured abundance [4] and collider searches [5]. In particular, we illustrate how the observation of displaced events at particle colliders could allow us to speculate about the early universe.
The main characters of our study are Feebly Interacting Massive Particles (FIMPs), new stable degrees of freedom beyond the Standard Model (SM) that never achieve thermal equilibrium through the expansion history of our universe. Nevertheless, they are viable DM candidates since decays and/or scatterings of primordial bath particles can produce FIMPs that free-stream subsequently. If FIMP interactions are renormalizable then their production is most efficient at temperatures around the mass of the heaviest particle involved in the process [6][7][8][9]. Consequently, the resulting relic density depends only on masses and couplings that we can measure in our laboratories and/or astrophysically and this "IR-dominated" production mechanism has been dubbed freeze-in in the literature [9]. Calculations for the relic density within this framework are reliable theoretically, but they come with an important caveat: the cosmological history is assumed to be standard. Contrarily, if the interactions proceed via non-renormalizable operators then the higher the temperature the more efficient production via scattering would be: the majority of FIMPs are produced around the time of inflationary reheating [10][11][12][13][14].
FIMPs are clearly a nightmare scenario for experimental searches. All we have left today are stable particles with expected rates at conventional DM searches -barring some exceptions, see e.g. [15] -typically too low to yield any signal. Here, we study a peculiar experimental manifestation of several FIMP models: new heavy particles, which participate in the FIMP production in the early universe, can be pair produced at particle colliders and then decay into DM with a macroscopic decay length. The resulting signature is the one of displaced decays with missing energy in the final state [16].
We investigate the scenario sketched in Figure 1: a FIMP X couples to a SM field A sm and a new beyond the SM (BSM) bath particle B via a cubic interaction. This will be the only DM coupling relevant phenomenologically. The FIMP is a SM gauge singlet and we take the cubic interaction to be very suppressed in order to prevent thermalization. In contrast, the bath particle B carries SM gauge quantum numbers and it is in thermal equilibrium with the primordial plasma at early times. Known examples for B can be supersymmetric partners (gluino, wino, squarks, etc.), but we take a bottom-up approach and we do not commit to any specific realization. Although we focus on a 3-body interaction, the generalization of our analysis to interactions involving more particles is straightforward.
Bath particles populate the early universe with X's via two different classes of processes. The cubic interaction itself is responsible for B decays, and it also mediates binary collisions once we account for interactions among primordial bath constituents. One example of such a collision is sketched in the right panel of Figure 1 where A sm and A sm are two SM particles that can differ from A sm . Production through decay is always IR dominated whereas the latter process can be UV dominated if the FIMP interaction is non-renormalizable. B < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > X < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > B < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " s e i W S e 4 F v r 4 1 3 n Y d E a h u O / x S 7 O g = " > A A A B / X i c d V D L S g N B E J y N r x h f U Y 9 e B o P g K c w E M c k t 6 E V v C Z g H J E u Y n X S S I b M P Z m a F s A Q / w K t + g j f x 6 r f 4 B f 6 G s 0 k E D V r Q U F R 1 0 9 3 l R V J o Q 8 i H k 1 l b 3 9 j c y m 7 n d n b 3 9 g / y h 0 c t H c a K Q 5 O H M l Q d j 2 m Q I o C m E U Z C J 1 L A f E 9 C 2 5 t c p 3 7 7 H p Q W Y X B n p h G 4 P h s F Y i g 4 M 1 Z q X P X z B V I k h F B K c U p o + Z J Y U q 1 W S r S C a W p Z F N A S 9 X 7 + s z c I e e x D Y L h k W n c p i Y y b M G U E l z D L 9 W I N E e M T N o K u p Q H z Q b v J / N A Z P r P K A A 9 D Z S s w e K 7 + n E i Y r / X U 9 2 y n z 8 x Y r 3 q p + K e n 7 S l j G K y s N 8 O K m 4 g g i g 0 E f L F 9 G E t s Q p x G g Q d C A T d y a g n j S t g H M B 8 z x b i x g e V s M t / v 4 / 9 J q 1 S k p E g b F 4 X a 7 T K j L D p B p + g c U V R G N X S D 6 q i J O A L 0 i J 7 Q s / P g v D i v z t u i N e M s Z 4 7 R L z j v X 8 v 9 l d w = < / l a t e x i t > X < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " / l f H M S a P 1 1 k 6 J x F w K N H t L f 5 i m p A = " > A A A B / X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o M e A F 7 0 l Y B 6 Q L G F 2 t j c Z M j u 7 z M w K Y Q l + g F f 9 B G / i 1 W / x C / w N J 8 k e N F r Q U F R 1 0 9 0 V p I J r 4 7 q f T m l t f W N z q 7 x d 2 d n d 2 z + o H h 5 1 d J I p h m 2 W i E T 1 A q p R c I l t w 4 3 A X q q Q x o H A b j C 5 m f v d B 1 S a J / L e T F P 0 Y z q S P O K M G i u 1 e s N q z a 2 7 C 5 C / x C t I D Q o 0 h 9 W v Q Z i w L E Z p m K B a 9 z 0 3 N X 5 O l e F M 4 K w y y D S m l E 3 o C P u W S h q j 9 v P F o T N y Z p W Q R I m y J Q 1 Z q D 8 n c h p r P Y 0 D 2 x l T M 9 a r 3 l z 8 1 9 P 2 l D G G K + t N d O 3 n X K a Z Q c m W 2 6 N M E J O Q e R Q k 5 A q Z E V N L K F P c P k D Y m C r K j A 2 s Y p P x V n P 4 S z o X d c + t e 6 3 L W u O u y K g M J 3 A K 5 + D B F T T g F p r Q B g Y I T / A M L 8 6 j 8 + q 8 O e / L 1 p J T z B z D L z g f 3 5 2 / l b w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 x v t Y S R h 9 w l U I j 7 j T K 2 N M x n w C P M = " < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 m c G J z / U w M 7 Y + 9 Q F 2 y 5 S m C I 0 + L s = " > A A A B 9 X i c d V D J S g N B E O 1 x j X G L e v T S G E R P w 2 Q h y y 3 i R W 8 R z A L J G H o 6 P U m T n o X u G j U M + Q 8 v H h T x 6 r 9 4 8 2 / s S U Z Q 0 Q c F j / e q q K r n h I I r s K w P Y 2 l 5 Z X V t P b O R 3 d z a 3 t n N 7 e 2 3 V R B J y l o 0 E I H s O k Q x w X 3 W A g 6 C d U P J i O c I 1 n E m 5 4 n f u W V S 8 c C / h m n I b I + M f O 5 y S k B L N 2 e D P r B 7 U D R W 3 u x k k M t b Z r 1 S r 5 Q s b J n W H A k p l 2 r V I i 6 k S h 6 l a A 5 y 7 / 1 h Q C O P + U A F U a p X s E K w Y y K B U 8 F m 2 X 6 k W E j o h I x Y T 1 O f e E z Z 8 f z q G T 7 W y h C 7 g d T l A 5 6 r 3 y d i 4 i k 1 9 R z d 6 R E Y q 9 9 e I v 7 l 9 S J w a 3 b M / T A C 5 t P F I j c S G A K c R I C H X Figure 1: Basic setup of this work: a FIMP X couples to the thermal bath via a cubic interaction with a SM particle A sm and a new BSM particle B. This sources X production via B decays (left panel) and binary collisions (one example in the right panel).
This setup opens up a promising way to probe FIMPs: B particles can be pair-produced at colliders via SM gauge interactions, and they can subsequently decay to FIMPs [17][18][19][20][21][22][23][24]. If the B decay length is too short, less than approximately 100 µm, the signal would be indistinguishable from conventional DM searches targeting missing energy. On the contrary, if it is larger than the detector size ( 10 m) we would observe a charged track associated to the B particle crossing the full detector; although this would be a spectacular manifestation of BSM physics, it has no manifest connection with the invisible universe.
The main focus of our study are signatures that we get when the decay length lies in between these two extremes. In such a range, we can observe charged particles decaying to SM and missing energy displaced from the primary interaction vertex.
What typical B decay lengths should we expect? If we neglect the mass of A sm , which is irrelevant to the phenomenology we study, the system in Figure 1 is described by three independent parameters: the masses of B and X and the decay length (or equivalently the lifetime) of B. For a given mass spectrum, the relic density constraint sets the B decay length. Unfortunately, the typical decay length is way too large to give rise to displaced events within the detector for B masses such that we can produce them at colliders [9]. This well known observation motivates a consideration. Let us assume that we observe displaced events and that we measure the B mass and decay length. In a quite optimistic scenario, let us even suppose that the decay happens via a 3-body interaction such as in Figure 1 and that we reconstruct the DM mass as well [25]. As we show in Section 2, this would point to a cosmological disaster: the relic density of X should be much larger than the observed DM abundance. Is there any way out? The answer is yes, and the reason is the caveat mentioned above: the relic density calculation assumes a standard cosmological history with freeze-in production during a radiation dominated (RD) era. However, the earliest time we can probe indirectly is the one of Big Bang Nucleosynthesis (BBN) when the temperature was T BBN MeV and the age of the universe was approximately one second. While it is reasonable to extrapolate the standard expansion history to earlier times and higher temperatures, this is an assumption not supported by any observation: we do not know what our universe looked like when it was younger than one second, see e.g.
Ref. [26]. Thus displaced missing energy at colliders could indirectly provide us with some hints about the early universe: DM production could have happened during a non-standard cosmology, see for instance Refs. [17-19, 24, 27-30] within the context of freeze-in.
The standard history cannot be extrapolated back in time indefinitely. At early times, our universe undergoes through a phase of very fast (exponential) expansion with energy density dominated by one (or more) scalar field known as the inflaton. Later on, the inflaton oscillates around the minimum of its potential and our universe experiences a phase of early matter domination (MD) during which inflaton decays to SM and BSM particles populate the primordial thermal bath. Eventually, when the age of the universe is around the inflaton lifetime, the energy budget is overtaken by radiation and we enter the RD epoch that lasts until the time of the (late) matter-radiation equality (T eq 1 eV). The highest temperature of such a RD epoch is conventionally called the reheating temperature, and we denote this quantity important for our study with the symbol T R . We investigate how such an early MD era, with reheating temperature T R below the B mass, alters the prediction for the B decay length at colliders. While we focus for concreteness on inflationary reheating, our results hold also for the last stages of a generic early MD epoch such as the one due to supersymmetric moduli [31][32][33][34][35][36][37] or evaporating primordial black holes [30].
We revisit the prediction for the B decay length during a RD epoch in Section 2, and we discuss how it changes for an early MD. For values of the B mass and decay length relevant to collider physics, and once we set the DM mass to the lowest value allowed by small scale structures constraints, we get an upper bound on the reheating temperature T R . In Section 3, we classify the simplest microscopic realizations. Finally, in Section 4, we first present the searches performed at the LHC experiments for long-lived particles in Section 4.1 and then illustrate in Section 4.2 how collider searches can bound the viable parameter space for these DM scenarios. Our analysis provides a map between unconventional DM collider signatures and properties of the early universe. Within our hypothesis of modified cosmological history due to late time reheating after inflation, this feature would be the reheating temperature T R . We defer technical details to appendices, and we conclude in Section 5.

Cosmological implications of the B decay length
In this section, we present numerical results for the FIMP relic density and we discuss how to understand them with the help of approximate semi-analytical solutions that capture the relevant physics. Technical details about the cosmological histories considered in our work and freeze-in production can be found in the Appendices A and B, respectively.

Freeze-in production via B decays in a RD universe
Every freeze-in calculation in a RD universe relies upon two assumptions. First, the reheating temperature is set to be much higher than the mass of any of the particles participating in the DM production processes. Second, DM production through inflaton decays or initiated purely gravitationally during inflation [38][39][40][41][42] is subdominant. We begin our analysis with the investigation of X freeze-in production via B decays in a RD universe and we also make these conventional assumptions. The rigorous way to track the X number density n X is by solving a Boltzmann equation, as we outline in Appendix B, and the output of such a procedure is illustrated in Figure 2. We show numerical solutions for a fixed B mass and three different values of its decay width Γ B . The variable on the vertical axis is the so called comoving number density, and it is defined as Y X ≡ n X /s where s is the entropy density of the radiation bath. The evolution "time variable" we employ is the dimensionless combination x ≡ m B /T where T is the temperature of the thermal bath. The comoving number density changes in a RD universe only if there are number changing processes. In other words, Y X is a convenient variable because it scales out the effect of the Hubble expansion.
Initially, the X comoving density grows with time and the slope of such a growth can be understood through a simple analytical estimate. The amount of X particles produced when the thermal bath had a temperature T during a Hubble doubling time, namely the time it takes for the universe to double its size, results in (2.1) Here, R B→AsmX (T ) is the (temperature dependent) rate of particle production and t(T ) is the age of the universe when the temperature is T . We estimate the former by multiplying the B rest frame decay width by an additional factor of m B /T to account for time dilatation due to the kinetic energy of the thermal bath. The latter is approximately the inverse Hubble parameter, t(T ) 1/H(T ), and for a RD universe we have H(T ) ∼ T 2 /M Pl . Thus the comoving number density scales with x (i.e., with the temperature) as follows (2.2) Figure 2 makes the IR domination manifest: the lower the temperature, the more efficient DM production becomes. This is only valid as long as B is present significantly in the bath since we hit the Maxwell-Boltzmann suppression once we reach temperatures  Figure 3: Inflaton (blue line) and radiation bath (orange line) energy densities as a function of the scale factor for inflationary reheating. In the top-right part of the plot we zoom in at early times and we show how the radiation bath energy density arises from inflaton decays.
around the B mass. Consistently with this physical argument, we observe how the three lines in Figure 2 reach a constant values for x larger than a few. In what follows, we refer to the temperature at which the FIMP production peaks around the B mass as T F I , or equivalently x F I ∼ O(1). The X comoving density freezes to the constant value A more rigorous calculation gives the expression in Eq. (B.15) of Appendix B. If we impose the further constraint that the FIMP X makes all the DM we get This relation connects the DM mass, the mass of the mother particle and its decay width. Keeping in mind that we work within the natural unit system ( = c = 1), this is the typical decay length for B decays at colliders after they have been pair-produced. Unless we consider a very light FIMP, the decay length is beyond the detector size (∼ 10 m) for B mass values relevant to collider physics. Displaced events seem somewhat unlikely within any FIMP framework. In fact, we will argue now how their observation at colliders (i.e., B decays inside the detector) could instead point to freeze-in during an early MD era.

A motivated alternative: early MD era during inflationary reheating
We provide in Appendix A the Boltzmann equations describing inflationary reheating. For a specific choice of inflationary parameters, Figure 3 shows how the inflaton and radiation energy densities depend on the scale factor a describing the expanding universe. Here, we use the semi-analytical solution given in Appendix A to describe the underlying physics. Initially, there is no radiation bath and the inflaton dominates the energy budget with a costant density E 4 I . At some point, which we take corresponding to a value of the scale factor a in , inflation ends and the inflaton begins damped oscillations around its minimum. The  Figure 4: Radiation bath temperature (left) and the entropy in a comoving volume (right) as a function of the scale factor. We set E I and T R to the same values chosen in Figure 3.
energy density stored in these oscillations red-shifts as non-relativistic matter, ρ M ∝ a −3 , and we enter an early MD epoch. Meanwhile, inflaton decays populate a newly formed radiation bath. The inflaton decays completely when the age of the universe is approximately its lifetime, and we recover the standard RD universe. The reheating temperature T R is the highest temperature ever achieved by the radiation bath during this RD epoch. It is instructive to investigate bath properties during reheating. We show in Figure 4 the evolution of two important quantities: the temperature (left panel) and the entropy in a comoving volume S(T ) = s(T )(a(T )/a in ) 3 (right panel). For values of the scale factor a < a in , the temperature and the entropy were both zero. As soon as inflaton oscillations began, the bath temperature had a quick growth reaching its maximum value T max when the scale factor did not even get twice as big, a max 1.5 a in . Up to order one factors, the maximum temperature as a function of the inflationary parameters reads T max M Pl Γ M E 2 I 1/4 with Γ M the inflaton decay width. Afterwards, the bath temperature decreases as follows where a R is the value of the scale factor when the bath temperature is T R (M Pl Γ M ) 1/2 . The decrease in Eq. (2.5) is slower with respect to the T ∝ a −1 of a RD universe, and this behavior persists as long as the inflaton field is around dominating the universe. The inflaton decays eventually, and we enter the RD phase for values of the scale factor a > a R . The entropy in a comoving volume grows during the reheating phase. We compute this quantity with respect to its value at the maximum temperature This entropy dump provides a dilution to the FIMP produced via freeze-in. Finally, it is useful to find an approximate expression for the Hubble parameter during reheating (2.7)

FIMP production during an early MD era
We revisit freeze-in via B decays when FIMP production peaks during an early MD era, i.e. when m B > T R . Our working assumption throughout this analysis is that the maximum temperature of the radiation bath attained during reheating, T max , is larger than m B . Thus the thermal bath generated during reheating begin its existence with a full relativistic abundance of B particles. Once we ensure that T max > m B the resulting FIMP relic density depends only on T R . The methodology to solve numerically the Boltzmann equation for freeze-in during an early MD era can be found in Appendix B. Here, we present numerical solutions and we understand the underlying physics thanks to analytical estimates. Figure 5 shows the evolution of the DM comoving number density for fixed values of the B mass and decay length, but different values of the reheating temperature. We can potentially reproduce the DM relic density upon choosing an appropriate value for m X . The purple line does not present any substantial difference with respect to the solutions in Figure 2. Indeed, the reheating temperature in this case is much larger than the B mass so freeze-in happens long after inflationary reheating is over and the cosmological background is a RD universe. An estimate for the relic Y ∞ X is still provided by Eq. (2.3). The red and the blue lines feature a quite different behavior. The DM comoving number density reaches a peak and decreases afterwards before settling down to its asymptotic value. Freeze-in production happens during an early MD epoch for both cases since T R < m B . We use the expression in Eq. (2.1) to estimate the amount of FIMPs produced during a Hubble doubling time when the universe had a temperature T Y prod making use of the Hubble rate given in Eq. (2.7). The above expression approximates well the blue and red curves of Figure 5 for temperatures larger than the B mass (x 1).
The initial growth scales as x 5 instead of the x 3 found for a standard cosmological history. The x 5 and x 3 dependencies are also well visible in the purple curve slope when going from FIMP production in an early MD to production in a standard RD era, namely in the transition between x < 10 −2 and 10 −2 < x < 1. Besides the different power-law behavior, the message about IR domination is confirmed and actually reinforced. The expression in Eq. (2.8) can only account for Y X (T ) at x 1 because B hits the Maxwell-Boltzmann suppression. Unlike the case of freeze-in during RD, we cannot estimate Y ∞ X by evaluating Y prod X (T ) at a temperature T F I of the order of the B mass because we still have to account for the subsequent dilution of Y X (T ) until the reheating time. The dilution, due to entropy released by inflaton decays, is well visible in the red and blue curves of Figure 5 for x > 1 up to x ∼ m B /T R . We quantify the dilution factor by the ratio of the entropy in a comoving volume between a production time at a given temperature T and the time when we can trust entropy conservation again (namely temperature below T R ) where we use the analytical estimate of Eq. (2.6). We define the ratio of the produced comoving FIMP density Y prod X (T ) of Eq. (2.8) and of the dilution factor D(T ) as which results from the combination of FIMP production in an early MD era at temperature T and subsequent entropy dilution between T and T R . We notice in Eq. (2.10) an extreme IR domination with a power-law scaling as 1/T 10 , and such an IR domination is the reason why the resulting relic density does not depend on T MAX . As for the RD case, FIMP production is more efficient at the lowest possible temperature compatible with the presence of B particles in the thermal bath. Thus the same argument used above would suggest that the relic FIMP comoving density is approximatelỹ Y X (T ) evaluated at a temperature T F I that is of the order of the B mass. However, Ref. [17] observed that FIMP production is still in action when we are in the tail of the Maxwell-Boltzmann distribution for B; this is just a consequence of the extreme IR domination observed in Eq. (2.10). The precise calculation of the relic density would require numerically solving the Boltzmann equation, and we provide this derivation in Appendix B. Here, we report the analytical estimate for the X comoving density Besides a quite large enhancement factor of 10 4 , freeze-in during an early MD era leads to a (T R /m B ) 7 suppression to the FIMP relic density that we observe today compared to the RD case of Eq. (2.3). The latter scaling is compatible with our estimate Y ∞ X ∼Ỹ X (m B ). Consistently with (2.11), we notice how the red and the blue lines lead to an asymptotic relic density that is suppressed with respect to the purple line by approximately three and seven orders of magnitude, respectively.

FIMP production from scatterings and UV sensitivity
As we illustrate in Figure 1, scatterings are an irreducible contribution to DM production within our framework. We complete our analysis by discussing this additional channel, and we show how its relative importance depends on whether the interaction between X and the thermal bath is renormalizable or not. Let d be the mass dimension of this operator. At high enough temperatures, larger than any mass of the particles involved in the collision, the scaling of the interaction rate with the temperature follows from dimensional analysis where Λ is the mass scale appearing in the FIMP interaction. For non-renormalizable interactions, d > 4, this is the scale suppressing the operator. The amount of FIMP produced is still given by Eq. (2.1), with production rate driven by scatterings as in Eq. (2.12). For a RD universe, the resulting comoving number density scales with the temperature as The value d = 4.5 divides two distinct regimes, and DM production is more efficient at high temperatures for d > 4.5. This case, where the interaction is non-renormalizable, is often referred to as UV freeze-in [43][44][45][46]. On the contrary, for renormalizable operators scatterings dominate over decays only if a larger number of scattering processes can contribute to DM production or when the mass splitting between B and X is small enough to suppress the decay rate (see e.g. Refs. [23,47,48]).
The abundance of DM produced through scatterings via non-renormalizable interactions cannot be computed in a RD universe. As it is manifest from Eq. (2.13), the production diverges as we go back in time (i.e., higher T ). However, the estimate in Eq. (2.1) is only valid as long as FIMPs are out of equilibrium, and eventually at high enough temperatures the rate would become so large that FIMPs thermalize. We need to know to what extent we can extrapolate the RD universe back in time in order to provide a precise calculation. In other words, we need a temperature UV cutoff corresponding to the highest temperature for the RD phase, and the reheating temperature provides such a UV cutoff.
The expression in Eq. (2.1) still accounts for the amount of FIMP produced in a Hubble doubling time when the temperature was T , and the scattering rate in Eq. (2.12) is still valid. However, as we have already seen in Sec. 2.3 for decays, there are two new ingredients: we need to use the Hubble parameter during reheating given in Eq. (2.7), and we need to account for dilution due to the entropy released by the inflaton decays. We find (2.14) DM production via scattering for renormalizable operators (i.e. d ≤ 4) is IR dominated for both RD and MD. On the contrary, in the range 5 ≤ d ≤ 8 DM production is dominated at small temperatures during reheating and at large temperatures during RD. In other words, during an early MD era it is mostly efficient either around T = T R or T ∼ m B , depending on the hierarchy between these two scales. For d > 8, DM production is UV dominated for both RD and MD. In the latter case, the resulting relic density is sensitive to T max and to the details of reheating. Interestingly, for d < 8 DM production via scattering depends only on T R and not on T max [17,43,49].
One concrete example of non-renormalizable interactions is for a Weyl fermion singlet FIMP χ. We consider the dimension 5 operator coupling χ to gluons and to a new fermion λ a in the adjoint representation of the SU (3) color gauge group where G a µν the gluon field strength. This is for instance the well known case of gluinogluon-goldstino coupling in supersymmetric theories [50]. Such an interaction induces DM production both via decays λ → G χ and scatterings. The rate for the former scales proportionally to the decay width of the mother particle, Γ λ m 3 λ /Λ 2 , and therefore production via decays is IR dominated even for non-renormalizable operators. The operator in Eq. (2.15) induces also DM production through t-channel scattering processes. Considering for example the t-channel exchange of a gluon, we have the matrix element where g s is the strong coupling, s and t are Mandelstam variables, and Π(T ) is a thermal mass that we insert to regulate the t-channel IR divergence.
Taking Π(T ) = g s T for illustration, Figure 6 presents numerical results for the operator in Eq. (2.15). Fixing the value of m B , we consider two different choices for the combination (Λ, T R ) in the two panels. For both cases, we take Λ to be much larger than the weak scale 1 and we choose numerical values such that we reproduce similar asymptotic DM comoving densities. The reheating temperature is always below the new physics scale, T R < Λ, in order to provide consistency with the contact interaction. As argued above, we do not need to specify the value of T max for d = 5 but it is implicit that we assume the hierarchy T max < Λ for the same reason. The left and right panels have reheating temperature larger and smaller than the B mass, respectively. We show the evolution of the production rates with respect to the Hubble parameter (upper part) and the DM comoving density (lower part) as a function of x = m B /T . We keep the two contributions, decays (orange lines) and scatterings (blue lines). Vertical dashed lines identify the two key temperatures when freeze-in production is mostly efficient around x F I ∼ O(1) (the purple lines are shown for x = 1 as a guide for the eye) and when inflationary reheating is over x = m B /T R (red lines). As expected, the scattering rate for the case T R > m B is mostly efficient around T R whereas the decay rate is maximized when the temperature is around m B . The relative height of the peaks is set by the reheating temperature: the higher T R , the higher the contributions from the scattering processes. Consequently, the comoving density features two plateaus corresponding to the two production channels. In contrast, both scattering and decay rates peak around T ≈ m B for the case T R < m B . For the specific interaction of Eq. the scattering rate is smaller than the decay rate around the peak and therefore scatterings provide a negligible contribution to the total DM abundance. Thus it does matter whether the operator mediating the decay is renormalizable or not because the resulting scattering production may have a relatively large importance and get UV dominated [10][11][12][13].

Learning about the early universe from displaced events
The main message from the analysis above is the one concerning the B decay length. The FIMP relic density is suppressed for freeze-in during an early MD, and consequently we need larger couplings than for freeze-in in RD. Thus if we want to account for the observed DM density we need a larger interaction rate, or in other words we need a shorter B decay length compared to Eq. (2.4): displaced events at colliders within the detectors is a rather typical prediction of freeze-in during an early MD epoch. Assuming that we can measure the masses from reconstructing the collider event in combination with information from measurements of the B production cross section -see e.g. Ref. [25] for new methods to read off the mass spectrum involved in displaced events -the observed decay length for B → A sm X will directly pinpoint the reheating temperature T R if we require that X makes all the observed DM. If we are less demanding and we require that X is just a subdominant component, it would at least provide an upper bound on T R .
Reconstructing the complete spectrum at collider is rather challenging so we should explore the less optimistic scenario where we cannot reconstruct the DM mass. FIMPs cannot be too light otherwise they would behave similarly to thermal warm DM (WDM) and erase small scales structures in contradiction with observations. Data from the Lyman-α forest put bounds on how warm these relics can be [51][52][53][54]. In particular, Ref. [54] put bounds on the WDM mass by implementing the velocity distribution into Boltzmann codes to get the linear transfer function and input the latter into dedicated hydrodynamical simulations. As a result, we have the bounds m wdm 3.5 keV (conservative) and m wdm 5.3 keV (stringent), corresponding to different assumptions on the thermal history of the universe. Refs. [45,[55][56][57][58][59] used the results of this analysis for WDM to obtain mass bounds on FIMPs, produced by decays or scatterings in a RD era, by comparing the associated linear transfer functions. 2 For the case of decays, assuming that the final state particles masses are negligible compared to initial state particle masses, Refs. [55,58,59] have obtained bounds in the range m X 7 − 16 keV. The strongest bound, obtained by Ref. [59] by employing the area criterium introduced in [61], results in m X 16 keV, and in the next sections we will refer to it as a "stringent" bound. As a guide for the eye, in the next sections, we will make use of the bound m X 10 keV as a "conservative" Lyman-α bound on FIMPs. It is worth keeping in mind that the freeze-in production via scatterings or decays during an early MD epoch might also affect these bound, see e.g. [45].
Even if we cannot reconstruct the DM mass, for a known B mass and decay length we can still determine a specific value of the reheating temperature T R once we set m X to the lowest value allowed by Lyman-α forest data. The DM mass cannot get any lower than that, and for such a smallest allowed value we need the smallest amount of dilution after freeze-in because ρ X ∝ m X . And the smallest amount of dilution corresponds to the largest value of T R ; this procedure also provides an upper bound on T R .
Finally, if we consider a non-renormalizable interaction DM production can be controlled by scatterings. If this is the case, the amount of DM grows with T R . If we fix the DM mass to the smallest value compatible with the Lyman-α bound, for a given m B and τ B we can derive an upper bound on T R , which is valid for any other (higher) DM mass. Thus, even when the dominant production channel would be freeze-in through scatterings, instead of decays, colliders could indirectly provide hints about cosmology.
To summarize, displaced events at colliders could indirectly provide constraints on the cosmological history of the early universe.

Minimal FIMP frameworks
The main goal of our analysis is to emphasize the interplay between early universe DM production and displaced events with missing energy at colliders. As argued in Section 2.5, a positive signal for this kind of events at the LHC would allow us to gain information about both the DM model parameter space and the cosmological history, even in the pessimistic case where the mass of the missing energy carrier can not be reconstructed.
We exploit this synergy considering a set of well defined scenarios giving rise to the cubic interaction shown in Figure 1. Restricting ourselves to BSM particles of spin 0 and 1/2, 3 and assuming DM to be a singlet under the SM gauge group, the resulting list is A sm Table 1: Classification of the simplest possible operators featuring a cubic interaction of the type of Figure 1. The DM particle is a SM singlet denoted by φ if it is a scalar, χ if it is a fermion. The bath field (Ψ B for fermions, Φ B for scalars) has therefore the same quantum numbers as the corresponding SM field A sm . See the text for further details.

Spin DM Spin B Interaction Label
compactly reported in Table 1. These "simplified models" are minimal in the sense that they feature the lowest dimensional operators (with d ≤ 5), containing one SM field and two extra fields of the dark sector, giving rise to the three-body interaction. Hence for each model only two fields need to be added to the SM field content. Furthermore, we impose an unbroken Z 2 symmetry under which the SM fields are even and the new fields are odd; this new parity, together with the fact that we always assume the DM field to be the lightest of the dark sector fields, ensures DM stability and restricts the three-field interactions to the form shown in Figure 1. In the first column of Table 1, we consider the finite set of options for the SM particle A sm : (i) a SM fermion field ψ sm , that is, a left-handed or right-handed lepton or quark; (ii) a SM gauge boson (photon, W , Z, gluon), following from an interaction involving the field strength The second and third columns contain all the corresponding possible choices for the spins of X and B, respectively. The gauge quantum numbers of B under the SM gauge group should be equal to the ones of A sm since the DM field is a gauge singlet. Notice how the case when A sm is the U (1) Y hypercharge gauge boson is of no phenomenological relevance at colliders because B would be a complete SM singlet for gauge invariance. In the fourth column, we display the Lagrangian term relevant for freeze-in DM production and for the decay of B into DM at colliders. We write these interactions denoting a fermionic DM particle (bath particle) as χ (Ψ B ) and a scalar DM particle (bath particle) as φ (Φ B ). In the last column, we give a label to each minimal model: the first letter indicates whether the bath particle B is a fermion (F) or a scalar (S), the second letter indicate the SM field involved, while the third letter indicate the nature of DM. It is worth noticing that all the interactions of Table 1 are renormalizable except for the ones of the models F F χ which are of dimension d = 5. The interaction term F F χ resembles the gluino-gluon-goldstino coupling or Wino-W -goldstino interactions in supersymmetry, cf. e.g. [50]. 4 Among the renormalizable interactions, three classes have d = 4. 4 One could in principle add analogous models featuring scalar DM and a scalar bath particle B in an adjoint representation of SU (2)L or SU (3)c with couplings to the field strengths given by dimension-6 operators of the kind ∂µΦB∂ν φF µν . However, as one can show by integrating by parts, such an interaction is proportional to the gauge boson masses, thus vanishing in the SU (3)c case and for the unbroken phase of SU (2)L (up to thermal corrections). We defer to future work a detailed study of this peculiar scenario. Some instances of models involving interactions of this kind have already been studied within the context of freeze-in in e.g. [48] for F ψ SM φ , in [8,24] for S ψ SM χ and in [21,22] for F Hχ . Finally, S Hφ is a dimension-three operator involving scalars only. This simple scheme encodes all possible cubic interactions involving one FIMP, another BSM field and one SM particle (restricting to spins 0 and 1/2 for the BSM sector). They amount to a total of 35 simplified models (counting all possible combinations of chirality and flavor of the SM fermions, and the three different choices for the gauge boson field strength F µν ). We will focus on a representative subset of these simplified models which can illustrate the general interplay between collider phenomenology and early universe cosmology.

Collider meets Cosmology
In recent years there has been an important effort of the LHC community to explore the sensitivity of the LHC to BSM signatures involving displaced decays and/or long-lived particles (together denoted as LLP searches) [16]. As argued in the introduction, LLP searches can test the nature of DM produced through freeze-in and they can also indirectly probe the cosmological parameters involved. In particular, the LHC LLP searches listed in Table 2 can directly probe the minimal frameworks presented in Section 3 as detailed in the summary Table 3. Within the context of an early MD era, we can obtain the DM relic density as a function of T R and confront the viable parameter space to the relevant collider searches. Below, we first underline the relevant features of the selected LLP searches in Section 4.1. In Section 4.2, we discuss how each of the theoretical models can be probed by experimental constraints on general grounds, and then illustrate the procedure for a representative set of minimal scenarios including leptophilic and topphilic as well as the nonrenormalizable case of the singlet-triplet model in Sections 4.2.1, 4.2.2 and 4.2.3 respectively. More details on the collider searches and our recastings can be found in Appendix C.

Collider searches for long-lived particles
We collect in Table 2 a selection of existing LLP searches performed by ATLAS and CMS employing data from collisions with center of mass (c.o.m.) energy of 13 TeV, and we assign a label to each of them. The selection is dictated by the attempt to have a full coverage over the possible SM final states with LLP signatures, and, when available, searches including missing transverse energy (MET) are preferred. In Appendix C.9 we comment about other possibly relevant LHC searches not considered in our study. Table 2 shows how ATLAS and CMS are spanning already a quite wide landscape of possible LLP-type signatures, involving leptons, gauge bosons, charged tracks, and also MET. As pictures can sometimes summarize information better than a long text we illustrate these searches in Figure 7 while we defer to Appendix C a brief description of each experimental analysis. Details on our recastings are also provided in length in the same appendix.
A large MET component is clearly a smoking gun signature of DM. Only two of these experimental analyses explicitly require missing energy in the final state. Nevertheless, the other searches are often sensitive to our simplified DM models due to a rather inclusive event selection. We remark that missing energy is not typically required in LLP searches  since the SM background is already very much suppressed by the displacement selection. However, the presence of a large MET in the targeted BSM signature could also be exploited for triggering the relevant events. Hence for the purpose of detecting DM models with LLP signatures, the MET requirement could be beneficial.

Early matter dominated era and collider signatures
We can confront the minimal simplified models of Table 1 to the series of LHC searches listed in Table 2 and discussed in Appendix C. The result of this comparison is illustrated in Table 3. We distinguish models involving the first and second generation fermions ( and q) from models involving the third generation (τ and t), because of their different collider signatures. We can exploit all the available collider searches focusing on models of the classes F W χ , F φ & S χ and F tφ & S tχ . These classes are highlighted in green. Another characteristic signature that can arise from the type of models under study is the one of the so-called "kinked tracks" (KT) included in the last column of Table 3. This signature corresponds to a long-lived charged mediator 5 B decaying inside the detector into DM and an electron or muon. The track due to the daughter particle appears as the continuation of the track of the mediator, however there is a kink at the point where the decay takes place. No dedicated and fully experimental analysis has been performed  Table 2 and Table 3 and described in Appendix C. The green area represents the tracker, the blue (cyan) region denoted by HCAL (ECAL) refers to the hadronic (electromagnetic) calorimeter, the gray boxes represent the muon spectrometer (MS). Dotted lines refer to DM, dashed lines refer to the LLP and they are darker for searches that require a charged track associated to it.
addressing this type of signature. We nevertheless include KT in Table 3 to highlight the range of physics targets that such a search would have; the reader interested in models of which kinked tracks would be a key signature can find an early discussion in Ref. [76]. Moreover, constraints on KT signatures can be obtained by re-interpreting the disappearing track searches along the lines of Ref. [77], as we will do in the rest of this section.
For our phenomenological discussion, we select within the model classes of Table 1 the following three that, we argue, represent an exhaustive set of simplified FIMP models: • S R χ , that features a scalar mediator with the quantum numbers of a right-handed lepton, coupling to DM through a renormalizable operator; • F t R φ , another example of renormalizable cubic interaction but involving stronglyinteracting fields, and the third generation (namely, a right-handed top); • F W χ , an explicit case of freeze-in controlled by a non-renormalizable operator, involving the SU (2) L field strength, hence a fermion mediator in a triplet representation.
Our selection includes renormalizable interactions (with first/second and third generation) and a non-renormalizable interaction. We perform a detailed analysis of the collider implication for early-universe cosmology for these models. As shown in Section 2, DM freeze-in production for renormalizable interactions will be IR dominated. In addition, for simplicity, we will assume that the DM only couples to one of the lepton/quark flavors at a time.
Notice that even for a mediator coupling to multiple families we do not expect any relevant constraint from flavor violation because of the smallness of the couplings considered, see e.g. [78] for a leptophilic scenario of such kind.

DV DJ DV
Label Table 3: Sensitivity of the LHC searches, shown in Table 2 and discussed in Appendix C, to our simplified models as labeled in Table 1. The fermion ψ sm of Table 1 can be either a light lepton ∈ {e, µ}, a tau lepton τ , a light quark q ∈ {u, d, s, c, b} or a top quark t.
Fermions can be either left-handed (L) or right-handed (R), a choice that has little impact on the final state but it can affect LHC production of the mediator Φ B /Ψ B . F W χ and F Gχ refer to models of the class F F χ where we consider the SU (2) L (coupling to electroweak gauge bosons) and the SU (3) c (coupling to gluons) field strength, respectively. The models we study in this work belong to the classes highlighted in light green.
The interplay between LHC and cosmology has been considered by Ref. [24] for models similar to our category F φ (thus with scalar instead of fermion DM) and F qφ (hence with first and second generation quarks involved, instead of the top). 6 Also, the model S t R χ was already studied in the context of freeze-in in [48], but without the special emphasis on the role of T R and on displaced events that we put here. In addition, R-hadron searches applied to a model F Gχ have been considered in [79] but in the context of freeze-out through coannihilations with feeble conversion processes (see also [80] in the case of S bχ ).

Leptophilic scenario, S R χ
We augment the SM field content by a new charged scalar Φ B with the quantum numbers of a right-handed lepton and a Majorana fermion singlet χ: this is the S R χ simplified model. Both extra fields are odd under an unbroken Z 2 symmetry. For simplicity and to ease the comparison of our results with a similar scalar DM model discussed in Ref. [24], we choose to couple the new fields to right-handed muons only, i.e. R = µ R . The Lagrangian reads where λ χ is a dimensionless Yukawa coupling. Within this model, there are two types of processes that contribute to the DM abundance via the freeze-in mechanism: (i) the decay of the mediator and (ii) the scattering processes converting the mediator into DM. 7 Notice that no symmetry prevents the new scalar Φ B to couple to the SM Higgs as In general, this quartic interaction will be radiatively generated at the level λ H 10 −2 , through e.g. a box loop involving electroweak gauge bosons. This interaction also contributes to the DM abundance through scattering processes involving H. However, we checked that its contribution is subleading, at most few percent for λ H ∼ O(1). We can thus safely neglect it, and in what follows we set the quartic coupling to λ H = 10 −2 .
Collider and cosmological constraints. The mediator Φ B can be pair produced at collider due to its electroweak (EW) gauge interactions. The production cross section at next to leading order (NLO) can be obtained by interpolating the tables available at [81]. We have checked that the NLO cross section is well reproduced by the leading order (LO) cross section obtained with MadGraph5_aMC@NLO [82] times a K-factor of approximatively 1.3. Depending on its lifetime, the mediator Φ B can lead to displaced leptons [68,69], disappearing/kinked tracks [65,66] and heavy stable charged particles [63,83] (see Table 3).
For long lifetimes, corresponding to a proper decay length exceeding the meter, we make use of CMS and ATLAS searches on HSCPs in Table 2, cf. Appendix C.1 for details of the recasting method. The corresponding regions excluded at 95% confidence level (CL) appear in red in Figures 8 and 9.
For intermediate lifetimes, when the decay of the charged mediators occur inside the tracker, one expects kinked tracks. Disappearing track searches can have some sensitivity to this type of signatures, that we estimate in the following. First, the charged track of the mother particle (the mediator) and the one of the daughter particle (the muon) should have a sufficient angular separation such that the track of the mother particle is interpreted as a disappearing track. In our analysis we conservatively require that the angular separation of the daughter track with respect to the mother track should be larger than ∆R > 0.1. Second, the DT searches include a lepton veto, and in this DM model the second part of the charged track is a non-soft muon, since the mass splitting between the DM and mediator is large ∆m m µ . If the muon is reconstructed, the event would be discarded. However, by the same requirement above on ∆R, the emitted muon does not generically point towards the primary vertex and hence we assume that it is not properly reconstructed and the lepton veto does not apply. 8 We refer to the Appendix C.2 for more details. In the figures of this section we will denote as "DT as KT" the sensitivity lines obtained by reinterpreting the DT search for a KT signature through the procedure explained above. For this purpose we have used the CMS search [67] since it is the one with the largest luminosity. The area that would be excluded by KT following our reinterpretation appears as a dark-yellow denselyhatched region in Figures 8 and 9. As can be seen, our reinterpretation demonstrates that a dedicated KT would nicely complement the other LLP searches for cτ φ ∼ 1 m. 7 For t-channel processes involving massive bosons and a lepton in the propagator, t-channel divergences should be regulated by thermal corrections. We expect that the regulated process contribute similarly to other scatterings, with a few percent correction in the abundance that we neglect in the leptophilic case. 8 We thank Steven Lowette for discussions on this point.
For shorter lifetimes, we employ the very recent ATLAS search for displaced leptons [70]. In the auxiliary material provided for this analysis, one can find the case of oppositely charged dimuons final state that can be applied directly to our model, leading to the exclusion region in Figures 8 and 9 colored with light green (tagged DL ATLAS).
For even shorter lifetime, we exploit the CMS displaced lepton analysis [68,69] to estimate the possible reach of an analogous search but targeting a different flavor configuration. The CMS analysis looks for events with a displaced eµ pair, while our final states consists of a pair of displaced muons since DM couples to muons only. The CMS results in [68,69] can not then be straightforwardly used for our leptophilic model. However, the CMS collaboration has provided the efficiency tables for displaced electrons and muons derived in these searches [84]. We optimistically apply the same efficiencies to a di-muon final states to obtain the sensitivity lines of an hypothetical experimental search, analogous to [68,69], but performed instead for a µ + µ − final state. The area that would be excluded by our reinterpretation of the CMS DL searches appear as a green loosely hatched region in Figures 8  and 9 (tag DL CMS). More details about our procedure are provided in Appendix C.3. 9 It is interesting to observe that displaced lepton searches of ATLAS and CMS covers slightly different regions of displacement.
In addition to collider constraints, the Lyman-α observations mentioned in Section 2.5 impose a conservative bound of m X 10 keV (and a stringent bound of m X 16 keV) and further limit the viable parameter space for DM. In the left panel of Figure 8, the light gray region is excluded by the conservative Lyman-α constraint. In the right panel of Figure 8, we also show our conservative m X 10 keV Lyman-α bound with a light gray area as a guide for the eye. In fact, a dedicated analysis of FI as non-cold DM arising from an early MD era is still missing in the literature. Finally, in some regions of the parameter space, DM may eventually achieve thermodynamical equilibrium and the freeze-in computations presented in Section 2 do not hold. We shade the corresponding parameter space regions in dark gray in Figures 8 and 9.
Results for fixed T R . The parameter space under investigation, including also the early MD cosmology, is four-dimensional and it is described by the mediator mass m φ and proper decay length cτ φ , the DM mass m DM and the reheating temperature T R . In Figure 8, we show slices that identify the cosmological history, or equivalently we fix the value of T R .
In the left panel of Figure 8, we consider a standard cosmological scenario with a rather high reheating temperature, T R = 10 5 GeV T F I . We show in the upper part with different colored lines, each one correspondent to a different value of the DM mass, the relation between the mediator mass and its decay length necessary to account for all the observed DM abundance. The colored areas tagged with the labels HSCP, KT, DL ATLAS and CMS are excluded or can be probed by the corresponding collider searches as discussed above. The only relevant searches for this case are those for charged tracks, such as HSCP, because mediator will typically cross the detector entirely. LLP searches targeting displaced vertices can essentially not test the viable DM space. Although the detection of a charged track associated to Φ B would be a spectacular signal of new physics, it would not let us . We identify regions excluded by HSCP searches from CMS and ATLAS (HSCP -red) and by displaced same-flavor lepton searches by ATLAS (DL ATLAS -uniform green). The possible reach from displaced sameflavor lepton by CMS is shown with a green loosely hatched region while our recasting of disappearing track searches in the case of a kinked track signature (DT as KT) is shown as a dark yellow densely hatched region. We denote with a light gray area regions excluded by the Lyman-α bound m DM 10 keV, see text for details. Finally, in the dark gray region DM thermalize in the early universe and therefore cannot be produced via freeze-in. learn much about the possible connection to DM or the physics of the early universe. A similar conclusion can be reached looking at the bottom part, where the colored curves are associated to different decay lengths and we visualize the (m φ , m DM ) plane.
In the right panel of Figure 8, we illustrate how this situation can be turned into a more promising one when considering freeze-in occurring during an early MD era with T R = 20 GeV. The upper part shows how displaced searches can probe DM with masses m DM 1 GeV. From the bottom panel, it appears that the decay length of the mediator ranges from the order of a millimeter to hundreds of meters while still accommodating the correct relic abundance when m DM > 10 keV and 100 GeV< m φ <600 GeV. Leptophilic Figure 9: Leptophilic DM with m χ = 1 GeV (left) and m χ = 10 keV (right). The upper panel shows constraints in the (m φ , cτ φ ) plane while the lower panel identifies the viable parameter space in the (m φ , T R ) plane. Solid lines reproduce the observed relic density, and the color code is the same as in Figure 8. The value of T R associated to a given (m φ , cτ φ ) in the right panel essentially provides an upper bound on T R for any dark matter mass as m χ 10 keV due to small scale structures constraints, see the text for details.
FIMPs can thus be probed at colliders through HSCP, KT and DL signatures. The HSCP searches probe charged mediator masses up to m φ ∼ 400 GeV and cτ φ 5 m. The recent DL ATLAS search has a very strong reach, excluding mediator masses up to ∼ 500 GeV for cτ φ ∼ 10 cm, and covering a large portion of the interesting displacements, i.e. 0.3 cm cτ φ 80 cm. For even shorter decay lengths, the DL CMS search, under the assumption that it can be performed for same-flavor leptons, as explained above, provide further sensitivity. For intermediated lifetimes, it is interesting to observe how a KT search (or any search sensitive to such signature) can close the gap between HSCP and displaced leptons searches.
Results for fixed m DM . Measuring the DM mass from a displaced event is quite challenging. However, we can still learn something about the early universe, and in particular extract an upper bound on T R , from a displaced event. We illustrate this point, and exemplify the discussion in Section 2.5, with a comparison in Figure 9 between a 10 keV DM scenario (right panel) and a case with heavier FIMP (left panel). In the upper parts, colored lines identify in the (cτ φ , m φ ) plane the values of T R reproducing the observed DM abundance. Likewise, in the lower parts they give the values of cτ φ (again reproducing the DM relic density) in the (T R , m φ ) plane. The colored areas show the regions probed by the collider searches as discussed above.
In the top-left panel of Figure 9, the case for m χ = 1 GeV, all continuous curves illustrate that freeze-in occurs in an early MD epoch, i.e. T R < T F I . In order to get a lower mediator lifetime for a given m φ we have to reduce T R , as also shown in the bottomleft panel. Reducing cτ φ corresponds to effectively enhancing the coupling constant λ χ in Eq. (4.1), which in turn increases the DM production rate and gives a larger freeze-in contribution. Therefore we need the early MD era to last longer in order to provide more dilution. The same panel also show that, for a given value of the mediator lifetime, larger m φ imply larger values of T R . Indeed, by increasing m φ we suppress Y ∞ X and we need a shorter duration of the early MD epoch to provide less dilution, see Eq. (2.11).
Similar conclusions can be drawn from the bottom-left panel of Figure 9, where we plot contours of fixed mediator lifetime accounting for the correct relic abundance in the plane (m φ , T R ). In the upper part of the plot, the contours become vertical since we recover freeze-in in a RD era T R T F I . Focusing on the bottom curves, an increase in the reheating temperature implies an increase in the decay length, and then a decrease in λ χ . A larger reheating temperature (with T R < T F I ) means less entropy dilution between T F I and T R , and in order to compensate for this effect there should be less DM produced during freeze-in, i.e. the decay width should decrease.
The picture changes when we consider the lowest value of the DM mass compatible with the Lyman-α constraint, m χ = 10 keV. Since the relic density is proportional to the DM mass, less dilution is needed to give displaced-vertex signatures at the LHC. Equivalently, the relic density can be reproduced for larger values of T R as it is well visible comparing the two panels of Figure 9. Furthermore, the colored curves in the top-right panel do not always satisfy T R < T F I , i.e. both cases of freeze-in occurring during an early MD era and a RD era are possible. Indeed, the contours for T R = 10 4 GeV and 10 3 GeV are superposed. In the bottom panels, the contours at fixed values of cτ φ become vertical for T R m φ . The contours in the right panel of Figure 9 are associated to lightest FIMP allowed by Lyman-α constraints, and they inform us on the maximal T R allowed compatible with freeze-in production (in particular, not larger than Ω χ h 2 = 0.12). For instance, suppose that a long-lived charged scalar decaying into a muon and MET is observed and its mass and decay length are measured to be m φ ≈ 500 GeV, cτ φ ≈ 1 cm. In such a case, even without a measurement of the DM mass, Figure 9 tells us that this particle can be a freezein mediator compatible with the observed DM relic density only if T R 100 GeV. Such information would have a profound impact on our understanding of the early universe (for instance inflation, baryogenesis, etc). Given that for heavier DM the upper bound on the reheating temperature only becomes more stringent, the discussion of the following models will focus on the m DM = 10 keV case only.
As a final remark, we comment on the possible impact of future experiments. First, assuming that a LLP search at the LHC remains background free, we can estimate the reach of HL-LHC (3000 fb −1 ) by (linearly) rescaling its sensitivity on the mediator production cross section with the luminosity. For instance, studying the cross section dependence on the mediator mass m φ , we can estimate that the DL ATLAS analysis will probe masses up to ≈ 900 GeV. This will considerably extend the explored parameter space of the model also to regions where lower T R upper bounds could be inferred (see Figure 9 top right).
On the other hand, we observe that future dedicated detectors targeting larger values of cτ φ (such as MATHUSLA [85] or CODEX-b [86]) will not constitute a very sensitive probe for this specific simplified model. First, the long-lived particle is electromagnetically charged and hence the regions of large cτ φ will be covered efficiently by HSCP searches at HL-LHC. Second, large cτ φ corresponds to very small couplings and hence to scenarios where dilution is not necessary in order to obtain the correct relic abundance. As a consequence, even in the event of a discovery, it will not be possible to set an absolute upper bound on T R in the region of large cτ φ , as it is visible from the "overlapping" contours in the top-right panel of Figure 9.
The above discussion changes significantly in the case of models with non-renormalizable interactions and UV sensitivity, as we will discuss at the end of Section 4.2.3, where future facilities dedicated to LLP signatures would be beneficial in exploring further the parameter space and the connection with T R .

Topphilic scenario, F t R φ
The second simplified model we study is the one with a real scalar DM particle φ and a vectorlike fermion mediator Ψ B with the quantum numbers of right-handed up-type quarks. We consider a "topphilic" scenario where BSM particles couple only to the top where λ φ is a dimensionless Yukawa coupling. As usual, both the DM field and the mediator are odd under an unbroken Z 2 symmetry, and we neglect possible extra contributions to DM production from a DM-Higgs portal coupling [87,88] so as to focus on the cubic interaction sketched in Figure 1. Decay and scattering processes can both produce DM in this scenario. Moreover, the higher (color) degrees of freedom and the larger gauge coupling are such that scattering processes are expected to play a more important role than in our previous leptophilic model, see also the discussion in [48]. As a consequence, for fixed values of m φ and T R , one will need a longer lifetime in order to account for the relic abundance in this topphilic scenario than in models of the class F φ (or S χ studied above). In Figure 10 we consider the topphilic scenario for a DM mass of 10 keV, a value that saturates the conservative Lyman-α constraint and allows to extract an upper bound on T R , as illustrated at the end of Section 4.2.1. In the left panel, colored lines identifies the value of T R needed to reproduce the relic density in the (m ψ , cτ ψ ) plane. Likewise, colored lines in the right panel show the values of cτ ψ needed in the (m ψ , T R ) plane. Comparing contours of Figure 10 with those of the right panel of Figure 9, we can see that longer lifetimes can be reached for fixed value of the mediator mass and T R in models F t R φ . Except from the LHC searches, this is the main difference between the leptophilic and topphilic scenarios. The dependence on the masses and temperature scaling are the same since for both cases all processes producing DM are renormalizable and thus IR-dominated.  Figure 10: Tophilic DM with m χ = 10 keV. Contours for fixed value of T R account for the DM relic abundance in the (m ψ , cτ ψ ) plane (left), and contours of fixed cτ ψ do the same in the (m ψ , T R ) plane (right). The value of T R associated to a given (m ψ , cτ ψ ) provides an upper bound on T R for any dark matter mass as m χ 10 keV due to small scale structures constraints. The blue shaded area is excluded by the DV+MET search, the red area by R-hadron searches, the blue dashed area by DV+µ search, the dark blue area by the DJ search, and the green area by the DL search.
Collider signatures and constraints. Looking at Table 3, one can infer that the topphilic model can lead to a rich variety of signatures at the LHC. The mediator decays into top whose decays can in turn give both fully hadronic and semi-leptonic signatures. In order to assess the most relevant searches, we first obtain the production cross section of the colored vectorlike fermion mediator Ψ B at the LHC multiplying the LO result of MadGraph5_aMC@NLO by a flat K-factor of 1.6. This K-factor reproduces well the ratio of LO to NNLO cross-section for squark pair production tabulated in [89] and we assume that the same K-factor can be applied to the case of colored vectorlike fermion mediator.
A long-lived Ψ B form an R-hadron (RH) before decaying [90]. If the decay occurs outside the detector, the R-hadron leaves a highly ionized track. The region excluded at 95% CL by RH searches according to our recasting is shown as a red area in Figure 10. To obtain these results, we employed the public code SModelS [20,91] and the one obtained from the "LLP recasting repository" [92], and applied the analysis done for a stop mediator in both CMS [63] and ATLAS [64] searches (see Section C.1) to our vectorlike top mediator as the spin of the produced particle is not expected to affect the sensitivity [90]. 10 For shorter decay lengths, the decay happens before the R-hadron leaves the detector with a top quark in the final state that itself decays dominantly into a bW pair. This in turn can decay hadronically to form jets, giving rise to displaced vertices (DV) or delayed jet (DJ) signatures. For details about these searches and our implementations we refer to, respectively, Appendix C.4 and C.5. Concerning the DV analysis, let us emphasize that the algorithm built to reconstruct the events in [71] combines displaced vertices that are located within 1 mm from each other. This is of particular relevance when the mediator couples to a top quark which dominantly decays to a b-quark and a W bosons. In the latter case, two displaced vertices possibly arise from the decay of the long-lived mediator since jets originating from a b-quark are displaced themselves (with a proper decay length of O(0.5) mm). Treating this properly is beyond the scope of this work. Therefore, we conservatively take only the tracks originating from the W into account for the reconstruction of the displaced vertex. The region excluded by the DV+MET search is shown with a light blue region while the region excluded by the DJ+MET is shown in dark blue. We see that the DJ+MET analysis can probe larger cτ ψ than DV+MET, see discussion in Appendix C.5. We also display with a dashed blue line the area probed by DV+µ (cf. Appendix C.6) which is similar to the one of DV searches. This captures the topologies where one of the top quark decays leptonically to a muon. The sensitivity is comparable to the DV searches since the suppression of the signal due to the leptonic branching ratio of the W is approximately compensated by the higher integrated luminosity employed in the analysis.
Alternatively, the two tops can both decay leptonically inside the detector leading to a displaced lepton pair. This signature can be covered by the DL searches [68][69][70]. We focus on the very recent ATLAS analysis [70] which employs the largest set of data at √ s = 13 TeV. In the Appendix C.3 we discuss the details of our recasting procedure, whose uncertainty is significant for some values of the LLP mass.
Nevertheless our analysis provides an indicative estimate of the LHC reach on the leptonic channel, which we display as a uniform green region in Figure 10. Even if the signal yield is limited by the leptonic branching ratio of the W boson, the search has a certain coverage of the parameter space, which is however superseded by the displaced jets + MET search for all ranges of cτ ψ .
LLP searches for the topphilic F t R φ model probe a larger parameter space than in the leptophilic S R χ of Section 4.2.1. The main reason is a larger LHC production cross section because the mediator is colored, and fermions display a higher production cross section than scalars. We see that RH searches probe the mediator masses up to m ψ ≈ 1.6 TeV for cτ ψ 30 cm. The DV searches cover a large range of parameter space corresponding to masses up to m ψ ≈ 1.9 TeV and cτ ψ 0.1 mm, partially overlapping with RH searches. Complementary constraints are obtained making use of DJ searches that cover 1.5 TeV m ψ 2.1 TeV and 30 cm cτ ψ 100 m. As in the case of the leptophilic model considered above, the plots in Figure 10 show how a measurement of an LLP mass and lifetime could be employed to bound the reheating temperature. Also for this model, the impact of HL-LHC searches will considerably extend the mass reach and hence the coverage of regions corresponding to low T R bounds. By naively scaling the sensitivity of the DV+MET search on the production cross section of the colored mediator with the luminosity, we estimate that HL-LHC could reach m ψ ≈ 2.8 TeV.

Non-renormalizable operators: the singlet-triplet model, F W χ
The two examples considered so far feature renormalizable FIMP interactions and therefore DM production is dominated by physics in the IR. We now discuss the F W χ model where interactions are non-renormalizable and scattering processes contributing to DM production are UV-dominated. Looking at Table 3, this model appears to be sensitive to a large set of existing LLP searches. The two BSM fields of the F W χ model are an extra singlet fermion DM and an SU (2) L triplet fermion mediator that we denote with and with masses m S and m T respectively. Both new fields are odd under an unbroken Z 2 symmetry and couple to the EW gauge bosons through the Lagrangian: where W a µν is the SU (2) L field strength and Λ is the scale of new dynamics responsible for generating the higher dimensional operator. More gauge invariant interactions exist at dimension five (for instance Hχ T Hχ S , see for details Appendix D). However, our goal is to explore interactions of the kind depicted in Figure 1, and we will assume that other possible dimension-five operators are suppressed, possibly by a higher mass scale, to the extent that they can be neglected for what concerns both DM production and LHC phenomenology. We only consider reheating scenarios with T max < Λ in order for the effective description in Eq. (4.5) to be valid during DM production. The singlet-triplet model has been also considered in Ref. [95] in the context of freeze-out DM production with compressed spectra.
We work in the m T m S regime and we neglect any Higgs portal coupling. As a result, singlet/triplet mixing will have minor influence and the lightest neutral fermion under the Z 2 symmetry, i.e. the DM candidate, is essentially the singlet fermion χ χ S with mass m DM m S . The mediator multiplet has a neutral component that is essentially the neutral component of the triplet Ψ 0 B χ 0 T with mass m T . The charged components, corresponding to the charged components of the triplet Ψ ± B = χ ± T , have a slightly higher mass m C = m T +∆m with a splitting induced by EW quantum corrections ∆m 160 MeV, see Appendix D for details. As a result, the mass hierarchy of the new fields in our analysis is always: DM production in the early universe. The processes reported in Table 5 of Appendix D, decays and scatterings, contribute to DM production. As discussed in Section 2.4, scattering mediated by non-renormalizable operators are UV dominated. The t-channel exchange of a photon in Ψ 0 B e − → χe − gives rise to a transition matrix ∝ (g/Λ) 2 (where g is the SU (2) L gauge coupling) that has a form similar to that shown in Eq. (2.16) given that all particles involved are much lighter than the mediator particle Ψ 0 B . 11 Similarly to Figure 6, we compare in the top panels of Figure 11 the production rate from scatterings (in blue) to the one of the decays (in orange) normalized by the Hubble rate as a function of x = m T /T . In the bottom panels, we plot the contributions of decay and scatterings (orange and blue dotted, respectively) to the total DM yield (continuous green). We fix m T = 1 TeV, and the prefactor 1/Λ reproduces the relic density constraint  Figure 11: Singlet-triplet model. Above: production rate/Hubble ratio for decay (blue) and scattering (orange). Below: DM yield. We consider T R > T F I (left) and T R < T F I (right). The dashed lines denote T = T R (red) and T = m T (purple) as a proxy for T F I . For these parameters, the observed DM abundance is reproduced for m DM = 10 keV.
for m DM = 10 keV; this gives a lifetime for the neutral mediator Ψ 0 B of cτ 0 = 3 m when T R T F I and cτ 0 = 0.3 mm when T R < T F I . The main contribution to DM production in the left panel comes from scatterings consistently with our discussion in Section 2. A larger reheating temperature implies an enhanced DM production, leading ultimately to a smaller DM-mediator coupling (i.e. a longer mediator lifetime) needed to account for the relic abundance. In contrast, in the right panel of Figure 11 we recover the yield dilution post freeze-in for x > x F I observed both for renormalizable and non-renormalizable operators in Figures 5 and 6. In this second case scattering processes play a subleading role.
We consider again m DM = 10 keV, and we show in Figure 12 contours of the values of T R giving the correct relic density in the (m T , cτ 0 ) plane (left panel) and cτ 0 contours satisfying the same constraint in the (m T , T R ) plane (right panel). In the left panel, we see how higher T R induces larger values of the proper mediator lifetime without reaching the saturation effect at large T R that was observed in the case of renormalizable operators (see the top-right panel of Figure 9 and Figure 10). Furthermore, cτ 0 contours in the right panel never become T R independent. These are direct implications of the UV-dominated scatterings at work when T R T F I for non-renormalizable operators, requiring smaller couplings (increasing cτ 0 ) when increasing T R .
Collider constraints. It is helpful to visualize the mass spectrum and the decay channels of the triplet components illustrated in Figure 13. The neutral heavy fermion Ψ 0 B can decay into DM emitting a γ or a Z boson through the dimension-five operator, and the decay width scaling as ∼ m 3 T /Λ 2 is typically macroscopic and in the parameter space relevant to freeze-in. The heavy charged fermion can either decay into the neutral heavy fermion plus a soft pion through the exchange of a W or directly into DM plus a W boson. The first decay mode is purely due to gauge interactions, while the second one is due to the dimension-five operator, hence it is controlled by Λ. Different Ψ ± B decay products give different signatures at the LHC. In the right panel of Figure 13, we illustrate how the different branching ratios (BR) of Ψ ± B decay products depend on the neutral mediator lifetime cτ 0 , which we take as a proxy for the scale Λ (as cτ 0 ∝ Λ 2 ). The Λ parameter indeed drives the relative importance of the Ψ ± B decays induced by the dimension-five operator relative to the gauge-induced decays. The decay length of the Ψ ± B is at most few cm, when the decay is dominated by the gauge-induced interactions. Shorter lifetimes can be obtained only when Λ is small such that the decay channels into DM dominate.
We distinguish in Figure 13 between the hadronic and leptonic subsequent decays of the Z and W bosons. For small cτ 0 (small Λ), Ψ ± B principally decays into χW ± . For larger cτ 0 (larger Λ), Ψ ± B mainly decay into a soft π ± and a heavy neutral fermion Ψ 0 B which in turn decays to χ + jets, leptons or γ. Let us also emphasize that, at the LHC, mediator pair production involves the production of charged fermions Ψ + B Ψ − B or the associated production of the Ψ ± B Ψ 0 B states, in all cases through s-channel electroweak bosons exchange. 12 We denote with colored areas in Figure 12 the regions excluded by LLP searches at LHC. For large values of the Ψ 0 B lifetime, the strongest bound come from searches for disappearing tracks (DT) [65,66] excluding the purple region up to m C m T 480 GeV.  Figure 13: Typical decay processes of the singlet-triplet model (left). Charged mediator branching ratio to the products listed in the legend as a function of the neutral mediator decay length cτ 0 as a proxy for the scale Λ for a mediator mass m T = 500 GeV (right).
This topology is relevant when producing a charged Ψ ± B decaying to π ± Ψ 0 B with the pion being too soft to be detected and Ψ 0 B too long-lived to decay inside the detector. If instead Ψ 0 B decays inside the detector, into DM plus a Z or a photon, DT searches can be no longer sensitive since extra hits could be recorded in other parts of the detector. In order to avoid this issue, we conservatively require that the decay of the neutral component Ψ 0 B occurs outside the tracker for the DT search to be applicable (see Appendix C.2 for details). This is the reason why DT searches appear to be sensitive only for cτ 0 1 m. In the plots of Figure 12, we show the DT exclusion region following from the ATLAS search. Indeed, ATLAS performs better than CMS for small decay lengths (see the discussion in Appendix C.2) which seems to be more relevant for the model under consideration.
Other signatures arise when the neutral component of the triplet decays inside the detector, or when the charged component decays directly to the singlet and a W boson. 13 If the Z or W boson decay hadronically, a displaced vertex (DV) or delayed jet (DJ) + MET signature could be observed. The heavy neutral Ψ 0 B is in some cases produced in the decay chain of charged Ψ ± B , which is itself long lived (with a decay length of at most few cm), together with a soft pion (see Figure 13). In this case, we assume that Ψ 0 B is emitted in the same direction as Ψ ± B , and hence we add the two displacements to define the total displacement of the resulting gauge boson. We consider all combinations of production and decay modes (with the corresponding branching ratios) to obtain the final states leading to displaced jets plus missing energy. Final states always contain two DM particles, consistently with Z 2 conservation, and two gauge bosons which can be either W or Z/γ (plus possibly soft pions). The relevant decay chains for these searches are the ones leading to jets (i.e. the hadronic decays of W or Z). The hadronic BR of Ψ ± B , as shown in 13 The charged component will only decay to the singlet when Λ (equivalently, cτ0) is small enough that This predominantly happen inside the detector, see the right panel of Fig. 13. Figure 13, depends on the neutral component decay length. The BR of Ψ 0 B into hadronic final states is instead independent on the value of its decay length. In Figure 12, the regions excluded by the DV [71] and DJ [72] searches are shown in light and dark blue, respectively. For moderate values of the lifetime (cτ 0 10 m) they put strong bounds on the mass of the triplet, excluding up to m T 1.3 TeV.
The other LLP signatures shown in Table 3 arise in this model. Leptonic decays of the Z and W boson lead to a rich variety of topologies, including displaced vertices + µ [73], displaced lepton vertices [74] and displaced leptons [68][69][70]. However, because of the comparatively small leptonic BR of the electroweak gauge boson, we expect these searches to be less constraining that the ones targeting hadronic final states plus missing energy. From Figure 12 we observe for instance that the displaced jets+muon search is once again slightly less constraining than the DV search, as already happened in the topphilic model. In Figure 12 we show also some representative sensitivity curves for purely leptonic LHC searches. In uniform green we display the reach of the recent DL ATLAS analysis [70]. Note that this search targets same-flavour (e ± e ∓ and µ ± µ ∓ ) as well as different flavours (e ± µ ∓ ) final states, which are all possible topologies for the singlet-triplet model (from leptonically decaying W and/or Z). By combining all these channels we obtain the sensitivity in Figure 12. The reach of this search cannot overcome the search focussing on displaced jets+MET, but it provides nevertheless a significant coverage of the parameter space of the model. For completeness we also show the reach of the the displaced lepton vertex search [74] (black region denoted as DLV in Figure 12), see Appendix C.7 for details. This search requires the two leptons to point to the same displaced vertex, hence effectively targets displaced Z for our model.
We have already discussed for the previous simplified models how the large luminosity of HL-LHC can further increase the mass reach of the LLP searches. Here we would like instead to emphasise the impact on this model of future detectors targeting particles with large (O(10 − 100) m) decay lengths, such as MATHUSLA [85] and CODEX-b [86]. First, in the singlet-triplet model the long-lived particle is neutral and can naturally reach the far detector, hence we expect that these facilities will significantly probe the large cτ region (a dedicated investigation as the one performed for the singlet-doublet model in Ref. [22] is left for future works). Second, the freeze-in DM production is UV sensitive. As a consequence, also for very small couplings (corresponding to large cτ ) the value of T R is set by the relic abundance constraint, as it can be seen by the contours in the left panel of Figure 12, and thus one can infer an absolute upper bound on the allowed T R . For instance, as we can see from the figure, observing a signal of a 1 TeV mediator with lifetime cτ 10 4 cm would imply, in the singlet-triplet model, an upper bound on T R of approximately T R 10 5 GeV. This would have relevant implications for well-motivated cosmological scenarios such as leptogenesis and for models of inflation.
Finally, comparing the results in Figure 12 with those for the other models in Figures 9  and 10, one can see how the relative sensitivity of different LPP searches differs in the three scenarios. This would provide an important handle to discriminate among possible models, in the (lucky) case of the observation of a pattern of signals from different LLP searches. The problem of determining the theoretical model underlying the LLP production is quite generic. For significant progress in this direction, see e.g. [96].

Summary and conclusions
Decades after its first convincing gravitational hints, we are still not aware about the microscopic particle identity of DM. This is a serious unresolved issue in particle physics and cosmology, and one of the most urgent questions to answer in fundamental interactions. When we look at a younger universe by detecting photons from distant objects, there is a limit on how much we can travel back in time. The universe was opaque to electromagnetic radiation before recombination, which happened when its age was approximately 380,000 years, and the CMB curtain hides whatever was behind. We can still look at earlier times by observing the BBN snapshot that features a radiation dominated universe only one second old. Looking at even earlier times is perhaps possible by observing spectral features in primordial gravitational waves, see e.g. Refs. [97][98][99][100][101][102], but at the moment we have no information about the pre-BBN universe. The detection of displaced events at colliders with missing energy in the final state can provide hints about both these open questions.
Focusing on the concrete framework illustrated in Figure 1, we have investigated the interplay between collider physics and early universe cosmology. We have considered DM particles with feeble couplings to the visible sector through the three body interaction in Figure 1 where B is a new BSM particle in thermal equilibrium at early times with the primordial bath. These scenarios are not accessible at direct and indirect DM detection experiments, and they are quite a nightmare for experimentalists. Nevertheless, they can give rise to spectacular signals at colliders where the bath particle B is pair produced and decays with macroscopic lifetime into DM. In the early universe, DM particles are produced through the freeze-in production mechanism (via decays or scattering). The decay length of the bath particle B, which one can measure from the displaced events at particle colliders, is directly linked to the DM relic abundance.
We have scrutinized this link with the aim to provide a deeper connection between the DM LLP signature and the cosmological history influencing the DM production. First, we have reviewed DM freeze-in production for a standard early universe with its energy budget dominated by a thermal bath of relativistic particles. The resulting prediction for the B lifetime is typically too long to lead to displaced signatures at colliders unless the DM is significantly light (few keV). Then we have discussed how a non-standard early universe history can affect the predictions for DM abundance. We have focused in particular on an early MD era resulting from a late inflationary reheating of the universe. Low values of the reheating temperature T R , as a consequence of the dilution of the DM relic density due to entropy release from inflaton decays, allow for higher values of the DM coupling to the SM and open the possibility to probe a larger part of the parameter space with LLP searches. We have improved previous studies in literature by going beyond the instantaneous reheating approximation and by considering both IR-and UV-dominated freeze-in. Remarkably, under certain generic assumptions, we have shown that DM LLP signatures can imply an upper bound on the reheating temperature of our Universe.
We have provided in Table 1 a systematic classification of simplified models giving rise to the cubic interaction displayed in Figure 1. After reviewing the relevant LLP searches within this context, we provided an in-depth study of the viable parameter space of three selected models that cover a wide spectrum of collider signatures (see Table 3). The first two involve renormalizable operators with DM preferably coupling to muons and top quarks respectively, while the last one features a non-renormalizable operator with preferred DM coupling to the EW gauge bosons. Processes driving freeze-in production are IR-dominated when renormalizable operators are involved and/or freeze-in mainly occurs via decays. Instead, when non-renormalizable operators are responsible for the three body interaction, freeze-in become UV-dominated for T R T F I with scatterings playing a major role. The interplay between DM relic abundance, affected by T R , and LLP searches is shown in Figures 9, 10 and 12. Our study highlights the model discriminating power of the LLP searches that have been performed by LHC collaborations. While the minimal models of Table 1 should not be regarded as "realistic" SM extensions, they still serve to illustrate how the underlying new physics could be identified in presence of a pattern of observed exotic signatures (cf. Table 3). They also represent useful building blocks that a complete theory may need to incorporate to account for signals for LLP at colliders.
A full reconstruction of the displaced events kinematics is highly unlikely. Notwithstanding the fact that we will probably be unable to measure the DM mass from these events, we have shown how displaced missing energy can still provide hints about the early universe. By setting the DM mass to the lowest value allowed by the Lyman-α forest bounds and imposing the relic density constraint, we have argued that one can still derive a generic upper bound on T R if the B mass and mediator lifetime could be measured at colliders by means of LLP searches. Colliders could thus indirectly provide a constrain on the early universe cosmology.
While we have focused on an early period of MD following inflationary reheating, our strategy can be applied to other scenarios with different early cosmological histories. One motivated example is the one of a "fast-expanding" universe [103] where the Hubble parameter has a stronger power-law temperature dependence, as it is the case for a kination phase [19,[104][105][106].
Our investigation highlights the importance of DM LLP signatures at colliders as a powerful strategy to search for feebly interacting DM. Indeed, beyond their promising discovery prospects in view of the HL-LHC, they can also represent indirect probes of the cosmological history of our Universe.

A Cosmological Histories
We summarize in this appendix the key features of the two cosmological histories analyzed in this paper: a standard radiation dominated universe and inflationary reheating. The expanding universe is described by the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric where physical length scales are proportional to the scale factor a(t). The Hubble parameter H(t) quantifies the expansion rate and it is set by the energy content via the Friedmann equation Here, we employ the reduced Planck mass M Pl = (8πG) −1/2 = 2.4 × 10 18 GeV. The species contributing to the total energy density ρ tot depend on the cosmological history under consideration. Consequently, the explicit functional form of the Hubble parameter in terms of the cosmic time t also depends on the cosmological history.

A.1 Radiation dominated universe
Once we consider a thermal bath of relativistic particles dominating the early universe, the energy density can be conveniently expressed as follows Here, T is the bath temperature and g * (T ) accounts for the effective number of relativistic degrees of freedom in the primordial plasma. All SM particles are relativistic for temperatures much higher than the weak scale, and this gives a contribution g sm * = 106.75. On the contrary, at temperature as low as the BBN epoch ( T BBN fews MeV, for a recent reassessment see [107]) we have g BBN * = 10.75. Thus the effective number of relativistic degrees of freedom changes approximately by a factor of 10 in the temperature range of interest. Another important quantity describing the thermal bath is its entropy density where g * s (T ) accounts for the effective number of degrees of freedom contributing to the entropy density. In our analysis we can safely take g * (T ) g * s (T ). The entropy in a comoving volume, S = sa 3 , is a conserved quantity for such a radiation dominated (RD) universe. This is not the case anymore once we consider inflationary reheating because inflaton decays dump entropy in the radiation bath. The conservation of entropy, dS/dt = 0, provides us with the useful relation It is important to remember that it is valid only if entropy is conserved.

A.2 Inflationary reheating
Once the inflationary phase is over, the inflaton field oscillates around its minimum and it decays to radiation, a process known as reheating. The intricate physics of reheating is accounted for by the compact set of coupled Boltzmann equations Here, ρ M is the inflaton energy density and Γ M its decay width. The factor of 3H in Eq. (A.5) shows how the energy density stored in the inflaton oscillations evolves as nonrelativistic matter (as opposed the factor of 4H in Eq. (A.6)). Our universe undergoes an early matter dominated (MD) epoch. Strictly speaking, Eq. (A.6) is only valid if g * = g * s = const; for the values mentioned above, the largest error we can potentially make is approximately 10%. As we discuss in our work for each case, the error due to this approximation is much smaller in our study and thus we safely employ it. We determine the cosmological background by solving numerically the Boltzmann equation system in our analysis. In this Appendix, we provide an approximate analytical solution that allows us to understand the behavior of the energy densities as a function of time. We start solving the system at the time t in when the inflationary phase is terminated and the initial radiation energy density is vanishing, ρ R (t in ) = 0; even if there was any radiation after the Big Bang, it would have been inflated away. On the contrary, the initial inflaton energy density is equal to ρ M (t in ) = E 4 I with E I the energy scale of inflation. If we consider our universe at times much earlier than the reheating lifetime, t Γ −1 M , the inflaton energy density just red-shifts as the inverse cubic power of the scale factor with a in is the value of the FLRW scale factor at the initial time t in . Meanwhile, the radiation bath begins being populated by inflaton decays and the approximate solution for the its energy density during this initial reheating phase reads [17] The radiation temperature is connected to its energy density via Eq. (A.2). Initially, there is not radiation at all and the initial growth is controlled by the second term in the square brackets of Eq. (A.8). The radiation energy density and its temperature reach a maximum for a value of the scale factor a max = (8/3) 2/5 a in . Neglecting the g * dependence on the temperature, the associated maximum temperature results in Afterwards, the decrease of the radiation energy density and temperature is controlled by the first term in the square brackets of Eq. (A.8). The radiation energy density features the peculiar dependence ρ R ∝ a −3/2 ; this has to be compared with the adiabatic expansion for a RD universe where the energy density decrease is faster, ρ R ∝ a −4 . Likewise, we notice that the temperature behaves as T ∝ a −3/8 when the universe expands and cools down; for a RD universe the temperature decreases with the scale factor is faster, T ∝ a −1 .
The approximation that lead us to Eq. (A.8) breaks down when most inflatons decays and the energy budget starts being dominated by the radiation bath: we recover the RD universe that serves as a background for BBN. The transition to a RD universe happens when the temperature of the thermal bath has a value T R satisfying the condition The so called reheating temperature reads The radiation temperature spans a potentially large range from T max down to T R , and they are connected via the relation T max (T R E I ) 1/2 .

B Boltzmann equation for DM production
The general form of the Boltzmann equation tracking the DM number density reads Besides the dilution due to the expansion, an effect controlled by the Hubble parameter H, the number density n X can change because of processes altering the number of X particles. This effect is accounted for by collision operator C α , and the index α runs over all possible processes. Here, we provide the collision operators relevant to our analysis. Although our work focuses on B decays to a 2-body final state, we present the collision operator for the more general case of n-body decays Here, B i denotes a generic particle belonging to the primordial thermal bath (either a SM or a BSM degree of freedom). The associated collision operator reads C B→XB 2 ...Bn = n eq B Γ B→XB 2 ...Bn where Γ B→XB 2 ...Bn is the partial decay width and K 1,2 are modified Bessel functions. The equilibrium number density for B can be obtained upon the phase space integration with g B the internal degrees of freedom (spin, color, etc.) of the bath particle. The last equality holds for a Maxwell-Boltzmann statistics, f eq B = exp[−E B /T ], and it leads to an exponential suppression in the non-relativistic regime (T m B ). For a generic binary collision where as usual B i denotes particles in thermal equilibrium at early times, the collision operator reads Here, σ is the Lorentz invariant scattering cross section as a function of the (squared of the) energy s in the center of mass frame. The lower integration extreme is and is set by a kinematical threshold, whereas the function λ is defines as follows

B.1 FIMP production during a RD era
Besides decays and scattering, the X number density gets diluted as n X ∝ a −3 because the universe is expanding. At the same time, the entropy in a comoving volume, S = sa 3 , is conserved for a RD universe. Thus in the absence of processes changing the number of X's, the comoving number density is a conserved quantity. Such a variable is convenient because it scales out the effect of the expansion, and it changes over time only if there are number changing processes.
It is convenient to introduce a dimensionless "time variable" x = M/T to solve the equation numerically. The overall mass scale M is purely conventional. For production via decays it is advantageous to set it to the decaying particle mass m B , whereas for production via scattering it is practical to set it to the heaviest mass involved in the process. Regardless of the specific choice for M , the Boltzmann equation reads This expression, written in terms of only dimensionless quantities, is the most general equation describing DM production during RD. We solve it numerically in our analysis. For a 2-body decay, with collision operator from Eq. (B.3), the equation reads We set the appropriate initial condition, Y X (x = 0) = 0, and we solve this differential equation by just performing an integral. The asymptotic value of the comoving number density, which we take as its value for very large x (i.e., small temperatures), results in . (B.14) The integral can be computed analytically if we set g * g * s const so that the number of relativistic degrees of freedom do not vary during the DM production era (see Refs. [23,108] for relaxing this assumption). We set them at the time when freeze-in is mostly efficient, namely when T m B (x 1), and we find For the case of FIMP production via scattering, with collision operator given in Eq. (B.6), we need to know how the cross section depends on the Mandelstam variable s.
Finally, we can evaluate the contribution to the Ω parameter defined as follows where we define the critical density ρ cr = 3H 2 0 M 2 Pl in terms of the current value of the Hubble parameter expressed in the conventional form H 0 = 100h km s −1 Mpc −1 .

B.2 FIMP production during an early MD era
The optimal numerical technique to solve the Boltzmann equation is different if freezein happens during an early MD era. The starting point is still Eq. (B.1), but using the comoving number density is rather inconvenient. Most of FIMPs are produced around the freeze-in temperature T F I > T R , and their number density dilutes afterwards with the scale factor as n X ∝ a −3 . At the same time, as long as the thermal bath temperature is larger than T R , inflaton decays dump entropy in the primordial bath and the entropy in a comoving volume sa 3 grows with time. As a consequence, the comoving number density Y χ decreases with time when the bath temperature spans between T F I and T R .
This suggests a new variable accounting for the number density of DM particles Once freeze-in productions runs its course, the quantity X remains constant until the present time. The Boltzmann equation tracking the FIMP number density becomes The cosmic time t is not an ideal evolution variable, the scale factor is quite convenient instead. We trade time derivative with derivative with respect to a by using the Friedmann equation in Eq. (A.1), and we find The inflaton and radiation energy density as a function of the scale factor can be determined from the procedure described in Appendix A. In particular, such a procedure provides us with a relation between the scale factor a and the bath temperature T which we can use to express the collision operators as a function of the scale factor. The asymptotic value of X , once we set X (a in ) = 0 as initial condition, results in At late enough times, once the thermal bath temperature drops below T R and the universe is dominated by radiation, we can trust entropy conservation again. In order to compute the FIMP contribution to the current energy density, we identify a late temperature T * < T R and we evaluate n X (T * ) and the entropy density s(T * ). The corresponding comoving number density: which is constant in the subsequent evolution of the universe, i.e. Y ∞ X = Y X (T * ). The energy density today thus reads ρ X (t 0 ) = m X Y ∞ X s 0 where s 0 is the current entropy density.  Figure 14: Comparison between our implementation of T R (continuous curve) and a cut in the time integration of the DM yield as in e.g. [23,24].

B.3 Comparison with previous studies
We conclude this appendix with a comparison between our analysis and previous studies about freeze-in via decays with subsequent dilution. As already explained above, we consider the full reheating process that happens after the inflationary phase of the early universe. In particular, we study the evolution of the thermal bath since its maximum temperature T max all the way down to T R . Our typical situation is when there is the hierarchy T max T F I T R , in such a way that freeze-in happens during an early MD era with the thermal bath filled of B particles with a relativistic abundance. Later on, until we reach the reheating temperature T R , DM particles produced via these decays get diluted.
An alternative methodology, followed e.g. by Refs. [23,24], is to neglect the extension in time of the early MD era. Within this approach, the two temperatures T max and T R coincide and this is achieved if the reheating process is instantaneous. On one hand, the amount of DM predicted suffers from a Maxwell-Boltzmann suppression of the B particles at temperature T R m B . On the other hand, the dilution due to inflaton decays after DM production is not present within this methodology.
We illustrate in Figure 14 how these two different approaches affect the estimation of the collider prospects. We plot with continuous curves the results obtained following our methodology within the context of our leptophilic scenario of Section 4.2.1. The dashed lines of the same color are obtained following instead the instantaneous reheating approximation. For renormalizable operators, this changes by a factor ∼ 2 at most the estimate of the mother particle lifetime. We expect thus that the analysis provided in Ref. [24] for scalar DM coupled to a vector-like heavy lepton mediator would give qualitatively equivalent results when following the methodology presented here.

C LHC searches for LLPs: description and recasting
In this appendix we briefly describe the LLP LHC searches tabulated in Table 2 and provide more details on our recasting procedure for each of them.

C.1 Heavy stable charged particle (HSCP) and R-hadrons (RH)
Description. When a charged mediator has a long enough lifetime to cross the detector completely (cτ B > O(10) m), it can leave a highly ionized track in the detector. This type of signature is usually referred to as Heavy Stable Charged Particle (HSCP). When the charged particle also carries QCD color, the long-lived particle is expected to hadronize, hence one refers to "R-hadron" (RH) searches. Searches of this type have been performed by both CMS [63,83] and ATLAS [64]. Both collaborations have publicly provided efficiency tables for the cases of slepton-like HSCP and gluino-or squark-like R-hadrons. The efficiency tables of the CMS search are gathered in the public code SModelS [91,109], which we use for recasting this search, see [47] for details. In order to reinterpret the ATLAS search, a public code is provided in the "LLP recasting repository" [92], a repository on GitHub holding various example codes to recast existing LLP searches. Notice that we use both codes to constrain R-hadrons originating from vectorlike quarks instead of squarks. In the latter case we assume that the differences in efficiency due to the difference in spin are minor.
As can be seen in Table 2, one would expect the ATLAS search to be slightly more constraining since the employed luminosity is about three times higher. However, there is an important difference between ATLAS and CMS searches, which unfolds itself when the search strategy is applied to look for charged particles/R-hadrons with a mass of the order of 100-200 GeV. The main difference between the two searches relies on how the signal regions (SR) are defined. Both make use of the reconstructed mass to define the SR, however the CMS SR are defined such that any heavy LLP/R-hadron will fit in at least one of the SR, while the ATLAS search only probes heavy LLPs/R-hadrons with a reconstructed mass larger than a certain threshold. 14 Hence, the CMS search is more sensitive in the lower mass region. Details on our recasting of these searches are provided in the following.
Recasting. Searches for highly ionized tracks originating from heavy stable charged particles (HSCP) and R-Hadrons (RH) have been both performed by CMS [63] and ATLAS [64]. As discussed above, the main differences between these two searches are the integrated luminosity and the sensitivity to the low mass range of the HSCPs/RHs.
In order to reinterpret the search performed by CMS, we employed a publicly available code named SModelS [91,109]. This code makes use of the available efficiency tables in order to find the upper limit on the cross section in any given model. However, the efficiency tables included in SModelS are built by assuming that the HSCPs/RHs are detector stable (i.e. decay way outside the detector). For HSCPs/RHs with an intermediate lifetime, it can happen that only a fraction of the produced particles traverse the whole detector. SModelS takes this effect into account by multiplying the efficiency of a stable HSCP/RH by the probability that the HSCP/RH with a specified lifetime traverses the detector completely (for more details on our treatment, see [47]).
In order the reinterpret the ATLAS search, we again made use of a publicly available code in the "LLP recasting repository" on GitHub [92]. This code does not rely on previously released efficiency tables, but rather makes use of Pythia 8 [110] to perform an eventby-event based analysis such that the dependence on the lifetime can be estimated with greater accuracy. The recasting strategy is based on the information published by the ATLAS collaboration, see [111]. For the leptophilic scenario described in Section 4.2.1, all necessary ingredients to extract the constraints are present in both codes. There is however one caveat for the reinterpretation of the RH searches for the topphilic scenario (Section 4.2.2). Both of these codes are able to reinterpret the existing searches only for squark-or gluino-like RHs. 15 Since we wanted to use these codes to constrain an RH originating from a vectorlike top quark, we made some small modification such that they would treat a vectorlike top RH as a stop-like RH. As the only difference between these two cases is the spin of the particle the RH is originating from, we do not expect that full implementation of a vectorlike top would have a large impact on our results.

C.2 Disappearing tracks (DT) and kinked tracks (KT)
Description. Disappearing tracks (DT) can arise in models with a charged mediator that has a proper decay length cτ B ∼ O(10) cm, i.e. smaller than the size of the inner tracking system of ATLAS and CMS, and its decay products cannot be reconstructed, either because they are neutral or they carry low momentum. Searches for disappearing tracks have been performed at both ATLAS [65] and CMS [66,67] with 13 TeV data.
In order to recast DT searches, we make use of the publicly available efficiency table provided by the CMS and ATLAS collaborations, see [47] for details. One important difference between the CMS and ATLAS search is the range of lifetimes they probe. The innermost tracking layers in the ATLAS detector have been upgraded and can focus more on very short tracks (tracklets). In practice ATLAS can probe shorter lifetimes τ ∼ O(1) ns while CMS is more sensitive to larger lifetimes with τ ∼ O(10) ns. Below, we reproduce the existing limits on the simplified model with a wino LSP as presented by CMS [67] and ATLAS [65], in order to validate our re-interpretations.
Note that these searches could also have some sensitivity for models exhibiting kinked track (KT) topologies. Indeed, let's assume that there is a displaced decay producing a stable SM charged particle plus DM. If the track of the charged decay product is not reconstructed, this signature would resemble a disappearing track. Some effort has been already done to verify the sensitivity of the disappearing track searches to a kinked track signature [18,24,77,112]. We remind here that DT searches impose a lepton veto and thus they cannot be employed straightforwardly for DM models involving a direct DM coupling to leptons and leading to displaced hard leptons in the final states. We discuss how we handle this issue in Section 4.2.1 to provide an estimate of the DT search reach on the KT signature which is present in that model.

Recasting.
Here we reinterpret the disappearing track (DT) searches by ATLAS [65] and CMS [67], performed for a supersymmetric model with long-lived winos. For that purpose, we apply the technique that some of us developed in [47] and applied to the DT searches 15 In SModelS, only efficiency tables for squark-and gluino-like RHs are present. Also Pythia 8 is only able to consider the formation RHs originating from squarks or gluinos.  Figure 15: Comparison between the experimental result (red) and our reinterpreted result (blue) for the latest DT searches performed by CMS [67] (left) and ATLAS [65] (right).
in [65,66]. We validated our technique for the DT searches used here [65,67] by reanalysing the case of the wino within our framework and comparing our results to the ones presented in the CMS and ATLAS papers. Figure 15 illustrates that we find a very good agreement.
The DT signature is a smoking-gun signature for the singlet triplet model studied in Section 4.2.3, which is closely related to the supersymmetric wino model, for which all DT searches are originally performed. The DT originates from the charged component of the triplet decaying into the neutral one in association with a soft pion that evades detection: For the charged track to fulfill the disappearance requirement stated in the CMS and ATLAS papers, no hits in the outer layers of the tracker have to be associated to the track. If the lifetime of the neutral component is short enough Ψ 0 B for it to decay inside the detector, hits in the outer layers of the tracker might occur. Therefore, we require in our analysis that Ψ 0 B to decay outside the tracker. In order to take this into account, we multiply the efficiency by the probability for Ψ 0 B to decay inside the tracker by taking l outer /γβ eff =0.36, as done in [20], where l outer is the outer radius of the inner tracker. Notice that the CMS search [67] considered here has been performed on a larger data set than the ATLAS search (L = 140 fb −1 compared to 36.1 fb −1 ). Naively, one would then expect CMS to set the strongest constraints on our models. As mentioned above though, the ATLAS detector has been optimized to study shorter disappearing tracks than the those studied by CMS, hence the CMS and ATLAS searches are complementary as one can see from our reinterpretation of Figure 15.
We have also estimated the efficiency of the CMS disappearing track search in constraining a kinked-track signature for our leptophilic model (see Section 4.2.1). For the DT search to be sensitive to a kinked track, we require that the angle between the mother and daughter particle tracks be large enough so as to avoid them to be identified as collinear (and hence being disentangled from one single track). We took the conservative limit on the angle ∆R ≡ (∆φ) 2 + (∆η) 2 > 0.1. Additionally, the CMS search imposes a cut on the amount of energy deposited in the calorimeters within a cone of ∆R < 0.5 about the disappearing track, which has to be less than 10 GeV. When the charged daughter particle is a muon, it will not leave a large energy deposit in the calorimeter, and hence, we do not need any extra requirements. In contrast, when the charged daughter particle is an electron or some hadronic final state, such an extra cut on the energy deposit in the calorimeter has to be imposed in order to estimate the sensitivity to the kinked track signature.

C.3 Displaced leptons (DL)
Description. In cases where the long-lived mediator decays to a lepton within the tracker, one would get displaced lepton signatures. Searches have been performed for events containing leptons with large impact parameters, i.e. displaced leptons, by CMS at both the 8 TeV [68] and the 13 TeV [69] run and by ATLAS at 13 TeV [70]. Both CMS searches are maximally sensitive to particles with cτ B ∼ O(1) cm while the ATLAS search is sensitive to slightly larger values of cτ B ∼ O(10) cm. In our study we consider both 8 TeV and 13 TeV CMS searches since the 13 TeV analysis has been performed on a small data set (corresponding to an integrated luminosity of 2.6 fb −1 compared to the 19.7 fb −1 employed by the 8 TeV search) such that the 8 TeV search will be in many cases more sensitive.
Let us emphasize that only displaced e ± µ ∓ pairs (coming from different decaying particles) are looked for by the CMS searches. In the cases where the DM couples to one lepton flavor only, this specific search does not a priori apply. In Section 4.2.1 we nevertheless used the available efficiency maps to provide a naive estimate of the possible LHC sensitivity to displaced µ + µ − pairs.
Instead, the DL ATLAS search defines three signal regions depending on the lepton flavours, one where a displaced e ± µ ∓ pair is observed and two others where same flavor leptons are observed (e ± e ∓ or µ ± µ ∓ ). Hence, it is straightforwardly applicable also to cases where the DM couples to one lepton flavor only. More details on the extracted sensitivities are provided below.
Recasting. The DL searches performed by CMS [68,69] look for events containing a displaced electron-muon pair that do not necessarily originate from the same vertex. The CMS collaboration has provided the model-independent efficiency tables for detecting a single displaced electron/muon as a function of its displacement and transverse momentum. By using these tables and applying the cuts reported in Table 4 one can estimate the modeldependent efficiency for the three signal regions (SR) that are defined for this search. These are based on the transverse impact parameter of both leptons (|d 0 |) and they are referred to as SR1 for 0.05 cm > |d 0 | > 0.02 cm, SR2 for 0.1cm > |d 0 | > 0.05 cm, and SR3 for |d 0 | > 0.1 cm. SR2 and SR3 are the most constraining as the SM background is virtually zero in these signal regions.
As mentioned above, we make use of the DL searches to estimate the possible LHC sensitivity to a displaced µ + µ − pair. Since the experimental collaboration has published the reconstruction efficiency for the electron and muon separately, we obtained a quite accurate estimate for the total event efficiency of a displaced µ + µ − pair. Notice that we have assumed zero background events in SR2 and SR3, as is the case for the e ± µ ∓ topology. Performing a full background estimation for the µ + µ − topology is beyond the scope of this work, but one should remind that a more realistic background estimation can possibly reduce the resulting sensitivity. 0.02 cm < |d 0 | < 2 cm 0.02 cm < |d 0 | < 10 cm ∆R > 0.5 ∆R > 0.5 Table 4: DL searches at CMS: Cuts on the transverse momentum p T , pseudorapidity η and transverse impact parameter d 0 of the electron/muon and the angle ∆R = ∆φ 2 + ∆η 2 between the electron and muon for the 8 and 13 TeV searches.
As mentioned above, the DL search performed by ATLAS [70] defines three different signal regions (SR-ee, SR-eµ and SR-µµ) and hence uses three different trigger strategies, namely single-photon, diphoton, and muon trigger. The photon triggers select events with an energy deposit in the electromagnetic calorimeter greater than 140 GeV (single-photon) and 50 GeV (diphoton) while the muon trigger select events with a signature in the muon spectrometer with p T >60 GeV and |η| <1.05. For each event, two signal leptons are defined as the leptons with the highest transverse momentum. These leptons must further have an impact parameter |d 0 | between 3 and 300 mm, a transverse momenta p T > 65 GeV and |η| < 2.5. Apart from these cuts, there are two extra event requirements. First, there must be a clear separation between both leptons, ∆R > 0.2 and none of the muons can be cosmic tagged. The latter consists in requiring that the two signal muons cannot be produced back to back, a case in which they are considered of cosmic origin and are hence removed. For events passing these selection cuts, the ATLAS collaboration also provides the model-independent reconstruction efficiencies for displaced electrons and muons [113] to be further applied. However, this information is only given for a benchmark with LLP mass of 400 GeV and proper lifetime of 1 ns. Nevertheless, we use for concreteness the same efficiency maps also for benchmarks with different LLP mass and lifetime, so we expect some discrepancy with the experimental results in certain regions of the parameter space.
In Figure 16 we show the validation of our recasting procedure for two representative models for which the full exclusion curves are reported in the auxiliary material of the AT-LAS paper, that is, production of right-handed smuon and right-handed selectron. We see that our recasting procedure reproduces very well the experimental exclusion curves around the benchmark point for which detailed informations on the electron and muon efficiency are provided by the ATLAS collaboration, as explained above. For masses significantly smaller or larger than 400 GeV the recasting is instead less precise. We nevertheless observe that, even if not quantitatively precise, our validation shows a qualitative agreement with the experimental results and hence we employ it in the main text to estimate the expected sensitivity for the topphilic and the singlet-triplet model. Righthanded Selectron exclusion Figure 16: Validation of the DL performed by ATLAS [70]. Left: right-handed smuon production. Right: right-handed selectron production.

C.4 Displaced vertices + MET (DV+MET)
Description. We call displaced vertices + MET the case where the long-lived mediator decays to DM and a color charged object before reaching the end of the tracker. A jet will be produced but is not reconstructed using standard jet-clustering algorithms. The events will be rather analysed by looking at the individual tracks of the jet originating from a displaced vertex, see [71] for details. In order to recast this search, we made use of the information provided by the ATLAS collaboration on the HEPdata page of the search [114]. We validated our recasting techniques by applying them on the model studied by the ATLAS collaboration and obtained very similar results, as discussed below.
Recasting. The ATLAS collaboration has released a search for events with displaced vertices in combination with missing transverse momentum in [71]. They look for events with at least one displaced vertex containing at least 5 tracks and missing transverse energy (MET) larger or equal to 200 GeV. The efficiency tables are publicly available on HEPData together with a procedure for recasting the search to other models [114]. Following this procedure, we reproduced the results of the ATLAS collaboration for long-lived gluinos decaying into quarks and a neutralino by doing an event-by-event analysis. Our results can be seen in Figure 17. A deviation between the ATLAS limit and our recasting appears for small values of the mass splitting between the gluino and the neutralino, ∆m/m 0.4. This is because we have omitted a cut, placed on the jets transverse momentum for 75% of the events, that only has a significant impact on the selection efficiency if the mass splitting is small. In our study, we do not explore compressed spectra and hence we can neglect this condition safely. We applied this recasting strategy to two of our models, the topphilic and singlet-triplet scenarios of Sections 4.2.2 and 4.2.3. Some special treatments were required for both of these models, as discussed in the following.
In the topphilic scenario, a displaced jet can arise from a long-lived R-hadron decaying into a top quark and the DM scalar. The top quark will predominantly decay to a bottom together with two other quarks (through an intermediate W boson). The quarks will promptly form a jet. The bottom quark on the other hand eventually decays to a slightly-  [71] by comparing the ATLAS results (red) with the results from our recasting (blue). Left: Upper limit on the production cross section vs. the gluino lifetime for a fixed gluino mass (1.4 TeV) and neutralino mass (100 GeV). Right: Upper limit on the production cross section vs the mass of the neutralino for a fixed gluino mass (1.4 TeV) and gluino lifetime (1 ns). Notice the discrepancy between the ATLAS limit and our result for a small gluino-neutralino mass splitting (∆m/m 0.4) which is due to the omission of one of the selection criteria (see the text for details).
displaced jet, with a displacement ∼ O(mm) due to the non-negligible b lifetime. For the DV+MET search, reconstructed vertices that are within a distance of 1 mm of each other are seen as one single vertex. Hence the tracks originating from the bottom and from the two other quarks will not necessarily be re-combined into one single vertex. Hence to be conservative, we discard any track originating from decay of the bottom quark in the recasting of the DV+ MET search (we also applied this conservative approach in the DV+µ search discussed in Appendix C.6).
In the singlet-triplet model, a displaced jet can originate from the decay of the charged and the neutral component of the triplet. The neutral components can be produced in two ways, through its gauge interactions or from the decay of the charged component. The decay of the charged component to the neutral one (together with soft pions) is slightly displaced, O(cm), and this displacement has to be taken into account in order to correctly interpret the efficiency tables for the singlet-triplet model. Due to the small mass splitting between the charged and the neutral component of the triplet, they will be more or less collinear. Hence, to estimate the total displacement of the displaced jet, we simply add the displacements of the charged and neutral components of the triplet.

C.5 Delayed jets + MET (DJ+MET)
Description. A delayed jet has been defined as a jet that is observed in the calorimeter of the detector at a later time than one would expect from a jet that is produced at the primary vertex. Such a time delay can be due to a heavy, long-lived mediator that slowly crosses the detector before decaying into a jet. This search relies on the timing capabilities of the calorimeter and, as a result, there is no need to reconstruct the displaced vertex.
Here we use the search for delayed jets + MET that has been performed by CMS in [72].
Notice that the delayed jet + MET search is able to probe longer lifetimes (cτ B ∼ 1 − 10 m) than the displaced vertices + MET search discussed above (more sensitive to decaying particles with cτ B ∼ O(1 − 10) cm) as the calorimeter lies further from the centre than the tracker. The DV+MET search instead is more sensitive to smaller lifetimes, and actually, the time delay will not be enough to distinguish between a displaced or a prompt jet. Hence, it is useful to consider both DV+MET and DJ+MET searches together given their complementarity.
Recasting. In order to probe larger values of the lifetime of long-lived particles decaying into jets, the CMS collaboration made use of the timing capabilities of the calorimeter to look for delayed jets [72]. They studied explicitly the case of long-lived gluinos decaying into a gluon and a gravitino. The efficiency tables for such a scenario are publicly available, starting from gluino masses of 1 TeV. Since these efficiency maps are blind to any jet activity arising from the primary vertex, we can use them directly for the scenarios studied here.
The gluino model studied in [72] always gives rise to two delayed jets while in the models listed in Table 1, it can happen that only one delayed jet arises. 16 In order to address this difference, we assume that the efficiency for an event with one delayed jet is small enough ( 1jet 1) so as to approximate the efficiency for an event with N delayed jets with With this approximation we can derive the efficiency for an event with one delayed jet, making use of the publicly available efficiency tables for long-lived gluinos (involving events with two delayed jets), and obtain the efficiencies for events with an arbitrary number of delayed jets. By weighting the derived efficiencies with the correct combination of branching ratios, we can obtain model-dependent efficiency tables for the models studied in this paper.
In the singlet-triplet model of Section 4.2.3, we have made extra assumptions in order to extract the corresponding efficiencies. As already mention above, e.g. in the discussion of DV+MET searches, the neutral component of the triplet can get an extra displacement from the decay of the charged component with cτ C ∼ O(cm). Since the CMS delayed jet search probes relatively long lifetimes of the neutral component (cτ 0 ∼ O(m) compared to cτ C ∼ O(cm) ), we assume that we can ignore the extra displacement arising from the decay of the charged triplet component. We have also assumed the charged component decaying into the neutral one and soft pions to one with a 100% branching ratio, i.e. BR(Ψ ± B → Ψ 0 B π ± ) = 1. Indeed, even though Ψ ± B can directly decay to DM, this decay width depends on the new physics scale Λ, while the decay to Ψ 0 B is driven by gauge interactions, i.e. independent of Λ. Since we probe large lifetimes with the DJ+ MET search, we are always in a regime where Λ is large such that the direct decay of the charged component to DM is suppressed compared to the decay to the heavy neutral triplet component. Stop mass (GeV ) Figure 18: Validation of the DV + µ search by ATLAS. The red curve is the ATLAS result, the blue curve is our reinterpretation using the available efficiency tables.

C.6 Displaced vertices + muon (DV+µ)
Description. Another type of displaced vertex search is performed by ATLAS looking for events containing displaced vertices together with a displaced muon track [73,115]. The search defines two orthogonal signal regions, one containing events triggered by missing energy, the other one containing events triggered by a high p T muon. In the context of this paper, where we focus on simplified DM models, the MET signal region is expected to be the most efficient. The displaced vertex is reconstructed in a very similar way as for the DV+MET search discussed in Appendix C.4, so we expect these two searches to constrain similar ranges of lifetimes. The main differences between the two searches are the requirement of a muon and the higher luminosity used in the DV+µ search.
Recasting. We use the search for displaced vertices and a displaced muon performed by ATLAS in [73] (with √ s = 13 TeV and L = 132 fb −1 ). It is particularly relevant for our simplified models as the long-lived particles are always produced in pairs and Z and W bosons potentially arising in their decays can partially decay to leptons. The ATLAS search defines two signal regions (SR): in the first SR (SR1) the event is triggered by a large amount of missing energy (E miss T > 180 GeV) while in SR2 the event is triggered by a track in the muon spectrometer (with p T > 62 GeV, |η| < 1.05). In addition, in SR2, a cut of E miss T < 180 GeV is imposed so as to define two orthogonal SRs. The experimental collaboration has provided trigger efficiencies, as well as muon and displaced vertex reconstruction efficiencies. We also make use of the cuts listed in the Tables 1 and 2 of [73] to reinterpret this search by doing an event-by-event analysis. We validated our approach by applying our recasting to theresults of this validation can be found in Figure 18.

C.7 Displaced dilepton vertex (DLV)
Description. This kind of searches targets long-lived neutral particles decaying into a lepton pair (µ + µ − , e + e − , or µ ± e ∓ ). In this work, we focus on simplified DM models where the connections between dark and visible sector is governed by a three-particle interaction (see Table 1). This interaction also governs the two body decay of the LLP in collider experiments. Since one of the two daughter particles of the LLP is always the DM, there will be only one SM particle involved in this decay. Therefore, a displaced lepton vertex will occur in the considered models only when the SM particle is a Z boson which decays to two leptons. The leptonic branching ratio of the Z boson is only about 10%, while the hadronic branching ratio is about 70%. Hence, we expect that the displaced lepton vertex search performed by ATLAS [74] will be generically less constraining than the DV+MET or DV+µ searches discussed in Sections C.4 and C.6 above. Nevertheless, since the ATLAS collaboration provides useful informations to reinterpret this search in a model independent manner [116], we recast this analysis as discussed in the following and we apply it when relevant for the considered models.
Recasting. Here we address the ATLAS search for a pair of oppositely charged leptons originating from the same displaced vertex [74]. The ATLAS collaboration has provided a document detailing how to reinterpret this search for any model containing long-lived particles decaying to oppositely charged leptons [116]. In order to calculate the event acceptance, one has to apply some simple kinematic cuts, among which the most relevant ones are that the invariant mass of the lepton pair has to be above 12 GeV and their displacement must be larger than 2 mm. The overall acceptance can be obtained by doing an event-by-event analysis.
To obtain the detection efficiency of a displaced lepton pair, the experimental collaboration provides two parameterisations. One for the R-parity violation (RPV) SUSY model and one for a Z toy model. There are two important differences between these models. First, in the RPV SUSY model, the LLP is a bino-like neutralino produced from the decay of a heavy squark. Due to the large mass splitting between the heavy squarks and the LLP, the LLP will have on average a much higher p T in the RPV model than in the Z model. As a result, the physical displacement in the detector will be much larger (due to the larger boost) for the same proper lifetime and the leptons arising in the decay will be also much more collimated in the RPV case than in the Z case. Another peculiarity of the RPV SUSY model is that the LLP decays to a lepton pair and a neutralino, such that the displacement vector pointing from the primary vertex to the secondary vertex is not parallel to the momentum vector of the lepton pair. This is also the case for the models studied in this work since the LLP will always decay to a DM + SM (including leptons). As a result, here we use the efficiency parametrisation provided for the RPV SUSY model. We validated the procedure discussed in [116] for the RPV SUSY model in Figure 19. Two different cases are illustrated, one where the LLP decays to either an electron pair or an electron-muon pair (left) and one where it decays to a muon pair or an electron-muon pair (right).

C.8 Delayed photons (Dγ)
Description. If the long-lived mediator decays to a photon plus DM, the expected LLPtype signature is displaced photons plus missing energy. The CMS analysis [75], looking for delayed photons, covers this topology. When a photon is emitted from a secondary vertex, it reaches the electromagnetic calorimeters (ECAL) at a later time compared to a photon produced at the primary vertex. The reason for the delay is twofold, firstly, there is m_sq=700GeV, m_n1=50GeV, Channel =emu/mumu Figure 19: Validation of the displaced lepton vertex search performed by ATLAS by comparing the upper limit on the cross section obtained by ATLAS (red) with the one from our recasting (blue) for the RPV SUSY model studied in [74].
a difference in trajectory, and secondly, the long-lived mediator traverses the detector at a lower speed. Another distinctive feature of a delayed photon is that it enters the ECAL at non-normal impact angles. These two features are used to distinguish a photon originating from a secondary vertex from one produced at the primary vertex.
The models that we investigate in Section 4.2 include this signature but with lower rate compared to other final states, so we do not expect that this could lead to relevant constraints on the parameter space of the models considered. Hence, we chose not to reinterpret this search. Moreover, little information for recasting is provided in [75], hindering a trustful reinterpretation.

C.9 Other searches
For completeness, we mention here other √ s = 13 TeV LHC searches that could constrain our models [117][118][119][120][121]. They all focus on displaced signatures involving jets, employing different search strategies and selection criteria. Since they do not target specifically models of DM with large missing transverse energy, we do not expect that they can test the parameter space of our models better than the DV+MET and DJ+MET analyses that we have re-interpreted in the main body of the paper. A possible interesting exception is given by the ATLAS analysis [118] which searches for displaced jets in the muon spectrometer and hence could possibly provide a slightly better coverage of the region with large cτ . At present such ATLAS search employs only 36 fb −1 , so we can expect only a moderate impact on our parameter space, but it will be interesting in the future to consider possible updated analyses employing a similar search strategy.
Another search that could constrain our models but does not focus on signatures involving jets is the ATLAS analysis [122] searching for tracks with a high ionization energy loss. This search typically targets models with heavy charged LLPs similarly to the searches for HSCPs and R-hadrons. However, this search does not require the track to traverse the detector completely, so that it targets values of the order of few meters. In principle, the search is sensitive to our leptophilic and topphilic models. However, we checked that in the leptophilic model, the mass range we considered (a few hundred GeV) is too low for the track of the mediator to be efficiently distinguished from a track of a high p T SM particle. In the topphilic model, on the other hand, the search is mostly sensitive to the same cτ range as the DJ search [72], but we expect the latter to have a higher reach in the mediator mass. This can be seen by comparing the maximal reach of the two searches for the model that they both employ to interpret their results (long-lived gluino production).

D Details about the singlet-triplet model
In Section 4.2.3, we have studied the singlet-triplet model, denoted by F W χ , which is an extension of the SM featuring a singlet and a triplet fermion, χ S and χ T . The most general Lagrangian one can write including interaction terms of dimension ≤ 5 reads where κ, κ ST , κ S , κ T and κ T are dimensionless coefficients and Λ is a common UV physics scale. In Section 4.2.3, we have effectively assumed that κ κ ST , κ S , κ T , κ T so as to focus on the cubic interaction of Figure 1. A complementary analysis focusing on the Higgs portal to DM for a large range of portal couplings can be found in [95].
If we neglect the terms that involve the Higgs field, i.e. setting κ ST = κ S = κ T = κ T = 0, no mixing occurs between χ S and χ 0 T and we can assume that the singlet field is the DM particle while the triplet plays the role of the freeze-in mediator, that is χ ≡ χ S , Ψ 0,± B ≡ χ 0,± T , as we did in Section 4.2.3. A mass splitting between Ψ 0 B and Ψ ± B always arises due to quantum corrections. As a result, the mass of the charged component can be written m C = m T + ∆m, with [123] where f (r) = r 2r 3 ln r − 2r + r 2 − 4(r 2 + 2) ln Hence the mass splitting between the neutral and the charged components of the triplet depends on the leading order mass m T . However, such a dependence is soft and for m T in the 100 GeV−1 TeV range, ∆m ≈ 160 MeV. As a consequence of this small mass splitting, the charged component of the triplet can possibly decay into the neutral component and a soft pion, with a decay width of [123] Γ(Ψ ± B → Ψ 0 B π ± ) = 2G 2 F f 2 π ∆m 3 π 1 − m 2 π ∆m 2 . (D.4)

Initial state
Final state Table 5: Decay and scattering processes, Initial state → Final state, contributing to the freeze-in production of DM for the singlet-triplet model.
The gauge interactions of the triplet and the cubic-interaction term proportional to κ 17 in Eq. (D.1) can be expanded as follows: This Lagrangian captures all possible interactions between the DM, the mediator and other SM particles which give rise to decay and scattering processes that could lead to freeze-in production of DM. The associated processes are listed in Table 5. 17 Notice that in the discussion of Section 4.2.3, we take κ = 1, that is, we reabsorb the coupling in the definition of the scale Λ.