Type III seesaw for neutrino masses in U(1)B−L model with multi-component dark matter

We propose a B − L gauged extension of the Standard Model where light neutrino masses arise from type III seesaw mechanism. Unlike the minimal B − L model with three right handed neutrinos having unit lepton number each, the model with three fermion triplets is however not anomaly free. We show that the leftover triangle anomalies can be cancelled by two neutral Dirac fermions having fractional B − L charges, both of which are naturally stable by virtue of a remnant ℤ2×ℤ2′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathbb{Z}}_2^{\prime } $$\end{document} symmetry, naturally leading to a two component dark matter scenario without any ad-hoc symmetries. We constrain the model from all relevant phenomenological constraints including dark matter properties. Light neutrino mass and collider prospects are also discussed briefly. Due to additional neutral gauge bosons, the fermion triplets in type III seesaw can have enhanced production cross section in collider experiment.


Introduction
Origin of light neutrino mass and dark matter (DM) in the Universe have been well known mysteries in particle physics for many decades [1]. While two mass squared differences and three mixing angles in the neutrino sector are very well measured [2], evidence suggests that DM gives rise to around 26% of the present Universe's energy density. In terms of density parameter Ω DM and h = Hubble Parameter/(100 km s −1 Mpc −1 ), the present DM abundance is conventionally reported as [3]: Ω DM h 2 = 0.120 ± 0.001 at 68% CL. Apart from this evidence from cosmology experiments like Planck, there have been plenty of astrophysical evidences accumulated for several decades. Among them, the galaxy cluster observations by Fritz Zwicky [4] in 1933, observations of galaxy rotation curves in 1970's by Rubin and collaborators [5], the observation of the bullet cluster by Chandra observatory [6] along with several galaxy survey experiments which map the distribution of such matter based on their gravitational lensing effects are noteworthy. While all these evidences are based on purely gravitational interactions of dark matter, there have been significant efforts in hunting for other types of interactions, possibly of weak interaction type, at several dark JHEP12(2019)109 matter direct detection experiments. Such interactions are typical of weakly interacting massive particle (WIMP) dark matter, the most widely studied beyond standard model (BSM) framework for accommodating dark matter. Here a DM candidate typically with electroweak (EW) scale mass and interaction rate similar to EW interactions can give rise to the correct DM relic abundance, a remarkable coincidence often referred to as the WIMP Miracle [7]. The very same interactions could also give rise to dark matter nucleon scattering at an observable rate. However, experiments like LUX [8], PandaX-II [9,10] and XENON1T [11,12] have so far produced only negative results putting stringent upper limits on DM interactions with the standard model (SM) particles like quarks. Next generation direct detection experiments like LZ [13], XENONnT [14], DARWIN [15] and PandaX-30T [16] are going to probe DM interaction rates very close to the rates at which coherent neutrino-nucleus elastic scattering takes place beyond which it will be extremely difficult to distinguish the neutrino or DM origin of nuclear recoil events. Similar null results have been reported by collider experiments like the large hadron collider (LHC) [17] as well as indirect detection experiments putting stricter upper limits on DM annihilation to standard model (SM) particles [18], specially the charged ones which can finally lead to excess of gamma rays for WIMP type DM. Though the null results reported by these experiments do not rule out all the parameter space for a single particle WIMP type DM, it may be hinting at a much richer dark sector. Though a single component DM is a very minimal and predictive scenario to begin with, a richer dark sector may in fact be natural given the complicated visible sector. There have been several proposals for multi-component WIMP dark matter during last few years, some of which can be found in [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. Going from single component to multi-component DM sector can significantly alter the direct as well as indirect detection rates for DM. Since direct and indirect detection (considering annihilations only, for stable DM) rates of DM are directly proportional to the DM density and DM density squared respectively, multi-component DM models can find larger allowed region of parameter space by appropriate tuning of their relative abundance. In addition to that, such multi-component DM scenarios often give rise to very interesting signatures at direct and indirect detection experiments [23,32,37,[52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69].
Similar to the BSM proposals of DM, there have been plenty of ways proposed so far, in order to accommodate light neutrino masses and mixing. While the SM can not generate light neutrino mass at renormalisable level, one can generate a tiny Majorana mass for the neutrinos from the SM Higgs field (H) through the non-renormalisable dimension five Weinberg operator [70] (LLHH)/Λ, L ≡ lepton doublet, Λ ≡ unknown cut-off scale. Dynamical ways to generate such an operator are classified as seesaw mechanism where one or more heavy fields are responsible for tiny light neutrino mass, resulting from the seesaw between electroweak scale and the scale of heavy fields. Popular seesaw models are categorised as type I seesaw [71][72][73][74], type II seesaw [75][76][77][78][79], type III seesaw [80] and so on. However, none of these models and their predictions have found any experimental verifications at ongoing experiments like the LHC [81]. In view of this, it is not really outlandish to consider frameworks where DM and neutrino finds common origin opening up the possibility to probe such models at different frontiers with more observables. One popular scenario that is built upon such objectives is the scotogenic framework, originally

JHEP12(2019)109
proposed by Ma [82], where particles odd under an unbroken Z 2 sector take part in radiative generation of light neutrino masses while the lightest Z 2 odd particle is the DM candidate. The same formalism was extended to two component DM scenario by the authors of [83]. While multiple copies of Z 2 symmetries remain ad-hoc in these models, one could find their UV completion with additional gauge symmetries. For example, in a recent work [47], two component fermion DM was proposed as a new anomaly free gauged B −L model, where B and L correspond to baryon and lepton numbers respectively. 1 The neutrino mass in this model originate from just two right handed neutrinos, often referred to as the littlest seesaw (LS) model [86], leading to a vanishing lightest neutrino mass. While TeV scale minimal type I or LS model has limited scope of being verified at experiments like the LHC due to the gauge singlet nature of additional fermions, presence of new gauge symmetries like U(1) B−L . Motivated by this, we consider another interesting possibility, within a gauged B − L model, where light neutrino masses arise from a type III seesaw scenario, along with the possibility of a two component fermion DM naturally arising as a possible way to make the model anomaly free. While implementing type I seesaw with three or two right handed neutrinos is a matter of choice (as both are allowed due to the unknown lightest neutrino mass) with the former not requiring any additional fermions (possible DM candidates) for anomaly cancellation, the implementation of type III seesaw in a gauged B − L model inevitably brings in new anomalies due to the non-trivial transformation of fermion triplet fields under the SU(2) L gauge symmetry of the SM. Out of several possible solutions to the anomaly cancellation conditions, we discuss one possible scenario which naturally leads to a two component fermion DM scenario. The U(1) B−L gauge boson apart from dictating the relic abundance of two component DM, also enhances the production cross section of heavy fermion triplets in proton proton collisions at the LHC due to the on-shell nature of mediating neutral gauge boson. We constrain the model parameters from all relevant constraints from DM sector as well as collider bounds and present the available parameter space. This paper is organised as follows. In section 2, we give an overview of gauged B − L model with different solutions to anomaly conditions including the one we choose to discuss in details in this work. In section 3, we discuss our model in details followed by section 4 where we mention different existing constraints on model parameters. In section 5 we discuss in details the coupled Boltzmann equations for two component dark matter in our model followed by section 6 where we mention briefly the details of dark matter direct detection. In section 7, we discuss our results related to dark matter relic abundance and constraints on the model parameters from all our requirements and constraints. We briefly mention some interesting LHC signatures of our model in section 8 and then finally conclude in section 9.
2 Gauged B − L symmetry Gauged B − L symmetric extension of the SM [87][88][89][90][91][92] is one of the most popular and very well motivated BSM frameworks that has been studied quite well for a long time. Since

JHEP12(2019)109
the SM already has an accidental and global B − L symmetry at renormalisable level, it is very straightforward as well as motivating to uplift it to a gauge symmetry and study the consequences. Compared to an arbitrary Abelian gauge extension of the SM, gauged B − L symmetry emerges as a very natural and minimal possibility as the corresponding charges of all the SM fields under this new symmetry are well known. However, a U(1) B−L gauge symmetry with only the SM fermions is not anomaly free. This is because the triangle anomalies for both U(1) 3 B−L and the mixed U(1) B−L − (gravity) 2 diagrams are non vanishing. These triangle anomalies for the SM fermion content are given as Remarkably, if three right handed neutrinos are added to the model, they contribute A New leading to vanishing total of triangle anomalies. This is the most natural and economical U(1) B−L model where the fermion sector has three right handed neutrinos apart from the usual SM fermions and it has been known for a long time. However, there exist non-minimal ways of constructing anomaly free versions of U(1) B−L model. For example, it has been known for a few years that three right handed neutrinos with exotic B − L charges 5, −4, −4 can also give rise to vanishing triangle anomalies [93]. It is clear to see how the anomaly cancels, as follows This model was also discussed recently in the context of neutrino mass [94,95] and DM [96][97][98][99][100] by several groups. Another solution to the anomaly cancellation conditions with irrational B − L charges of new fermions was proposed by the authors of [101] where both DM and neutrino mass can have a common origin through radiative linear seesaw. Very recently, another anomaly free U(1) B−L framework was proposed where the additional right handed fermions possess more exotic B − L charges namely, −4/3, −1/3, −2/3, −2/3 [102]. The triangle anomalies get cancelled as follows.
These four chiral fermions constitute two Dirac fermion mass eigenstates, the lighter of which becomes the DM candidate having either thermal [102] or non-thermal origins [103].
The light neutrino mass in this model had its origin from a variant of type II seesaw JHEP12(2019)109 mechanism and hence remained disconnected to the anomaly cancellation conditions. In a follow up work by the authors of [104], these fermions with fractional charges were also responsible for generating light neutrino masses at one loop level. One can have even more exotic right handed fermions with B − L charges −17/3, 6, −10/3 so that the triangle anomalies cancel as In the recent work on U(1) B−L gauge symmetry with two component DM [47], the authors considered two right handed neutrinos with B − L number -1 each so that the model still remains anomalous. The remaining anomalies were cancelled by four chiral fermions with fractional B − L charges leading to two Dirac fermion mass eigenstates both of which are stable and hence DM candidates.
To implement type III seesaw in U(1) B−L model, we consider two copies (n Σ = 2) of fermion triplets Σ Ri (i = 1, 2) into the model having quantum numbers (1, 3, 0, −1) under SU(3) c , SU(2) L , U(1) Y and U(1) B−L gauge groups respectively. Note that, two copies are enough to satisfy the light neutrino data, as in LS models. The non-vanishing anomalies are Note that the first anomaly is arising only due to the SU(2) L transformation of newly introduced fermions and was absent in the minimal B − L extension of SM. Obviously, the first anomaly can be cancelled only by introducing an additional field which has nontrivial transformation under both SU(2) L and U(1) B−L . We introduce such additional fields with a goal to keep our setup minimal and connected to the origin of light neutrino mass. Introducing a quintuplet Ψ(1, 4, 0, n 1 ), the first anomaly becomes which can vanish if n 1 = 2n Σ /5 = 4/5. The other anomalies can be cancelled by introducing additional fields which do not contribute to the first anomaly and hence SU(2) L singlets. The remaining anomalies are These remaining anomalies can be cancelled by introducing three SU(2) L singlet chiral fermions , N 2L 1, 1, 0, 2 5 , N 3L 1, 1, 0, 6 5 . (2.7) -5 -

JHEP12(2019)109
This can be seen as follows Since the quintuplet and the other singlet fermions have no role to play in generating light neutrino mass, we do not pursue this possibility further. We can have three fermion triplets: two of them having B − L charge −1 and the third having exotic charge n 1 . We will check if the third fermion can have any possible role in generating light neutrino masses. In such a case, there arises a possibility to get vanishing [SU(2) L ] 2 U(1) B−L anomaly. In this case, which can vanish if n 1 = n Σ = 2. The other anomalies can be cancelled by introducing additional fields which do not contribute to the first anomaly and hence SU(2) L singlets. The remaining anomalies are We now consider different possible solutions to these anomalies one by one. This can be seen as follows However, this solution is not very motivating owing to the existence of singlet fermions having B − L charge −1, which will give type I seesaw contribution to light neutrino masses, already discussed by several earlier works.
Solution 3: the remaining anomalies can also be cancelled by the following fermions with fractional B − L charges: (2.10) In order to have non-zero masses for all new fermions and sticking to minimal scalar contents having integer B − L charges, we find that there can be one stable dark matter candidate, in terms of one of the singlet fermions. Since single component dark matter in such models have already been discussed in several works, we do not pursue this possibility further.

JHEP12(2019)109
Solution 4: the most interesting possibility is the solution which can also be recast as We can construct two Dirac fermions from these four chiral ones, just by introducing two singlet scalars φ 1 , φ 2 having B − L charges 1, 4 respectively. The corresponding mass terms will be Another scalar singlet having B − L charge 2 is also required in order to give mass to the fermion triplets Σ 1,2 . The third fermion triplet acquires mass and also mixes with the other two fermion triplets by virtue of its couplings to the scalars φ 1 , φ 2 . Although the third fermion triplet is not stable due to its mixing with the first two, the two Dirac fermions constructed above can be separately stable and hence give rise to multi-component dark matter. Also, the third fermion triplet, through its mixing with the first two, can contribute to the light neutrino mass matrix as we discuss in details in upcoming sections. Due to these interesting possibilities, we pursue this scenario in our work.

The model
In this section we have discussed our model elaborately. As mentioned earlier, in this work our prime motivation is to have a multicomponent dark matter scenario where both dark matter candidates are stable by virtue of a single symmetry group. Additionally, we want to address neutrino mass generation as well. Keeping these two things in our mind, we have extended the SM in all three sectors namely the gauge sector, the fermionic sector and the scalar sector. In the gauge sector, we have demanded an additional local U(1) B−L gauge invariance where B and L are denoting baryon and lepton numbers respectively of a particular field. This introduces anomalies (both axial vector and gauge-gravitational anomalies) in the theory which can only be evaded by the inclusion of additional fermionic degrees of freedom. This has elaborately been discussed in the previous section. We have seen that for the case of three additional SU(2) L triplet fermions Σ iR , i = 1 to 3, two of which with B − L charge −1 are necessary to generate neutrino masses via Type-III seesaw mechanism and the rest having B − L charge 2 is solely required to cancel [SU(2) L ] 2 U(1) B−L anomaly. We further need four SM gauge singlet chiral fermions with fractional B − L charges to cancel both U(1) 3 B−L and (gravity) 2 U(1) B−L anomalies. Moreover, at least three scalar fields φ i (i = 1 to 3) are also necessary to give masses to all new fermions in the broken phase of the U(1) B−L symmetry. We have properly adjusted B − L charges of these φ i s such that only the Dirac mass terms among these singlet fermions are possible and more importantly the Dirac mass matrix is diagonal. This results in two physical Dirac fermions out of these four chiral fermions which are simultaneously stable and thus both can be viable dark matter candidates. In tables 1 and 2, we have listed all fermions as well as JHEP12(2019)109 (1, 1, 0, − 14 5 ) Table 1. Fermionic fields of the present Model including the SM fermions.
The Lagrangian of our present model invariant under the full symmetry group is given by Here, L SM denotes the SM Lagrangian involving quarks, gluons, charged leptons, left handed neutrinos and electroweak gauge bosons. The second term is the kinetic term of From table 2, we have already understood that our model has a very rich scalar sector and the gauge invariant interactions among the scalar fields are described by L scalar which contains following terms,

JHEP12(2019)109
where covariant derivatives for the Higgs doublet H and singlet scalars φ i s are defined as The quantity D H µ is the usual covariant derivative of the SM Higgs doublet with g and g are gauge couplings of SU(2) L and U(1) Y respectively and the corresponding gauge bosons are denoted by W a µ (a = 1, 3) and B µ . The covariant derivative of H does not include B − L gauge boson Z BL as the corresponding gauge charge of H is zero. On the other hand, being a SM gauge singlet, covariant derivative of φ i only contains Z B−L with n φ i is the respective B − L charge of φ i and g BL is the new gauge coupling. After breaking of both B − L symmetry and electroweak symmetry by the VEVs of H and φ i s, the doublet and all three singlets are given by where v and u i s (i = 1, 2, 3) are VEVs of H and φ i s respectively. For calculational simplicity we have assumed all three VEVs of singlet scalars are equal i.e. u 1 = u 2 = u 3 = u. Therefore, substituting eq. (3.4) in eq. (3.3) we have a 4 × 4 mixing matrix for the real scalar fields in the basis 1 √ 2 (h s 1 s 2 s 3 ) T which has the following form, The physical scalars (h, s 1 , s 2 , s 3 ) are obtained by diagonalising the above real symmetric mass matrix and they are related by an orthogonal transformation to the unphysical scalars (i.e. the basis state before diagonalisation) as where O is a 4 × 4 orthogonal matrix which makes M rs diagonal i.e. O T M rs O ⇒ a diagonal matrix containing all the masses of four physical scalars as diagonal elements.
In appendix A, we have expressed all the elements of O matrix in terms of four mixing angles (assuming mixing among the three singlet scalars are identical). Additionally, for simplicity we also set β = 0, 2 and ζ = δ in eq. (3.3).
2 β is the coupling of quartic interaction among all four singlet scalars which has less significant impact in dark matter phenomenology compared to other cubic scalar interaction terms like δ φ1φ1φ † 3 and ζ φ3φ3φ † 2 .

JHEP12(2019)109
On the other hand, in our model we have four pseudo scalars as well, which are z, A 1 , A 2 and A 3 (see eq. (3.4)). Out of these four pseudo scalars, z does not mixes with others and becomes the Goldstone boson corresponding to the SM Z boson after EWSB. However, the pseudo scalars of complex singlet φ i s mix among each other when B − L symmetry is broken by the VEVs of φ i s. Therefore, unlike the case of real scalars, here we have a 3 × 3 mixing matrix in the basis 1 After diagonalising, this matrix we get two physical pseudo scalars A 2 , A 3 and a massless Goldstone boson A 1 corresponding to the extra neutral gauge boson Z BL . This can easily be checked as the pseudo scalar mass matrix has a null determinant. The eigenvalues of the pseudo scalar mass matrix for ζ = δ, β = 0 are 3 where ζ must be less than zero to ensure m A 2 and m A 3 are real. Further, the physical CPodd scalars (A i s) and the Goldstone boson A 1 are related to the unphysical basis states Moreover, the neutral gauge boson Z BL becomes massive after the breaking of U(1) B−L and the corresponding mass term of Z B−L is given by Now, let us concentrate on the fermionic sector of the present model. Here, in addition to the usual SM fermions, we have four gauge singlet fermions and three SU(2) L triplet fermions. All these new fermions have appropriate B − L charges. In eq. (3.1), L f is the Lagrangian for these newly added fermionic fields and it is composed of two parts as

JHEP12(2019)109
where, the Lagrangian for the singlet fields are given below In the above Y i (i = 1, 2) are the dimensionless Yukawa couplings and the covariant derivative is defined as is the corresponding B − L charge of N κL(R) which is listed in table 1. As mentioned earlier, due to special choice of B − L charges of scalar fields φ 1 , φ 2 and φ 3 , the Yukawa interactions in eq. (3.12) are exactly diagonal in the basis ξ 1 = N 1L + N 1R and ξ 2 = N 2L + N 2R . In this basis above Lagrangian can be rewritten as, which can be further simplified as where P L,R = 1 ± γ 5 2 , left and right chiral projection operators. Besides, in the above we have assumed the Yukawa couplings Y i s are real. Therefore, from the above Lagrangian one can easily notice that both ξ 1 and ξ 2 are decoupled from each other and thus can be stable simultaneously. Hence, they naturally form a two component dark matter system stabilises by the B − L symmetry only. 4 On the other hand, SU(2) L × U(1) B−L invariant 5 Lagrangian for the triplet fields are given by, where, we consider fermion triplets Σ kR and its CP conjugate Σ kR c in 2×2 representation as where C is the charge conjugation matrix. The covariant derivatives used in the kinetic terms of Σ kR and Σ c kR can be defined as where n k Σ is the B − L charge of Σ kR and there is a sign flip as the B − L charges of Σ kR and its CP conjugate are equal but opposite in sign. Now, we define ψ 0 Majorana fermion and a Dirac fermion ψ − k = Σ − kR + Σ + kR c . Following [105], we have written the triplet Lagrangian in terms of ψ 0 k and ψ − k as θ W = tan −1 g g is the weak mixing angle (Weinberg angle). The last two terms of the above Lagrangian introduce off-diagonal elements in the mass matrices of both ψ 0 k and ψ − k (k runs from 1 to 3) when s 1 gets a nonzero VEV. As a result, one needs to diagonalise both the mass matrices using bi-unitary transformations in order to get the physical fermionic states. However, in this work we have not considered this. We have worked in a limit when two triplet fermions having B − L charge −1 will play crucial role in generating observed neutrino masses and mixings via Type-III seesaw mechanism. The Yukawa Lagrangian involving the leptons is given by The light neutrino mass matrix is generated from the type III seesaw mechanism as where the Dirac neutrino mass matrix M D and neutral triplet fermion mass matrix M R are given by Diagonalisation of the light neutrino mass matrix using the above forms of M D , M R gives one vanishing mass eigenvalue. This is same as the prediction of the recent work [47] as well as the littlest seesaw model [86] mentioned earlier. If we simplify our neutral fermion triplet mass matrix M R by incorporating the smallness on off-diagonal Yukawa couplings and equality of singlet VEVs, we can approximate M R to be a diagonal matrix leading to a relatively simple light neutrino mass matrix given as This simplified light neutrino mass matrix also leads to a vanishing lightest neutrino mass while the three mixing angles can be satisfied by suitable tuning of the Yukawa couplings. Since for TeV scale triplet fermions, the Dirac Yukawa couplings y αβ Σ have to be fine tuned at the level of ≤ 10 −4 , they do not impact the dark matter analysis discussed in this work. Hence, we do not take such couplings into our subsequent discussions.

Constraints on the model parameters
Before going into the detailed calculation of DM relic abundance and relevant parameter scan, we note the existing theoretical as well as experimental constraints on the model parameters.
In order to keep the scalar potential (given in eq. (3.3) within the square brackets) bounded from below, the quartic couplings must satisfy the following inequalities: To prevent perturbative breakdown of the model, all quartic, Yukawa and gauge couplings should obey the following limits at any energy scale: Experimental limits from LEP II constrains such new gauge sector by putting a lower bound on the ratio of new gauge boson mass to the new gauge coupling M Z /g ≥ 7 TeV [106,107]. The corresponding bounds from the LHC experiment have become stronger than this by now. Search for high mass dilepton resonances have put strong bounds on mass of such gauge boson coupling to first two generations of leptons with couplings similar to electroweak ones. The latest bounds from the ATLAS experiment [108,109] and the CMS experiment [110] at the LHC rule out such gauge boson masses below 4-5 TeV from analysis of 13 TeV data. Such bounds get weaker, if the corresponding gauge couplings are weaker [108] than the electroweak gauge couplings. Also, if the Z gauge boson couples only to the third generation of leptons, all such collider bounds become much weaker, as explored in the context of DM and collider searches in a recent work [111].
Similarly, the additional scalars in the model also face stringent constraints which typically arise due to their mixing with the SM Higgs boson which enable them to couple with the SM particles. The bound on such scalar mixing angles would come from both theoretical and experimental constraints [112,113]. In case of scalar singlet extension of SM, the strongest bound on scalar-SM Higgs mixing angle (θ 1j , j = 2, 3, 4) comes form W boson mass correction [114] at NLO for 250 GeV M s i 850 GeV as (0.2 sin θ 1j 0.3) where M s i is the mass of other physical Higgs. Whereas, for M s i > 850 GeV, the bounds from the requirement of perturbativity and unitarity of the theory turn dominant which gives sin θ 1j 0.2. For lower values i.e. M s i < 250 GeV, the LHC and LEP direct search [115,116] and measured Higgs signal strength [116] restrict the mixing angle sin θ 1j dominantly ( 0.25). The bounds from the measured value of EW precision parameter are mild for M s i < 1 TeV. While these constraints restrict the singlet scalar mixing with SM Higgs denoted by (θ 1j , j = 2, 3, 4), the other three angles (θ 23 , θ 24 , θ 34 ) remain unconstrained. We choose our benchmark values of singlet scalar masses and their mixing with SM Higgs boson in such a way that these constraints are automatically satisfied.

The Boltzmann equations for two component DM
In this work, as we already know that we are dealing with two stable dark matter candidates ξ 1 and ξ 2 . To find the present number densities of dark matter candidates we need to solve JHEP12(2019)109 two Boltzmann equations one for each candidate. The collision term of each Boltzmann equation contains all possible number changing interactions of that particular dark matter candidate allowed by the symmetries. In the present model, there are two types of number changing interactions for a dark matter candidate. First one is the pair annihilation where ξ i andξ i annihilate in a pair into all possible final states (X) except a pair of other dark matter candidate ξ jξj (j = i). These processes reduce the number of ξ i andξ i by one unit (assuming there is no asymmetry in the number densities of dark matter and its antiparticle). The other type of number changing process is ξ iξi → ξ jξj (i = j). This is actually the conversion process where one type of dark matter converts into another. It increases the number of lighter dark matter candidate by two unit while reducing the number of heaver one by the same amount. This conversion process acts as a coupling between the individual Boltzmann equations of ξ 1 and ξ 2 . Let n 2 = n ξ 2 + nξ 2 and n 1 = n ξ 1 + nξ 1 are the total number densities of two dark matter candidates respectively. Assuming there is no asymmetry in number densities of ξ i andξ i , the two coupled Boltzmann equations in terms of n 2 and n 1 are given below [23,117,118], where, n eq i is the equilibrium number density of dark matter species i. The second term in the left hand side involving the Hubble parameter H represents dilution of number density due to the expansion of the Universe. The extra half factors in the collision term are due to non-self-conjugate (Dirac fermion) nature of our both dark matter candidates ξ 1 and ξ 2 [119]. Here we are more interested to study the evolution of number densities of two component WIMP system due to all possible interactions which change particle numbers of either ξ 1 or ξ 2 or both. Hence, in stead of actual number density n i , it is convenient to describe the Boltzmann equation for a species in terms of comoving number density Y i = n i /s with s representing the entropy density of the Universe. The quantity Y i is an useful one as it absorbs the effect of the expansion on n i . In terms of comoving number densities, the above two coupled Boltzmann equations can be written as,

JHEP12(2019)109
where, G is the Gravitational constant and x ξ i = m ξ i T , is a dimensionless variable with T being the temperature of the Universe. The quantity g is expressed as, where, h eff (T ) and g eff (T ) are the effective degrees of freedom related to the entropy and energy densities of radiation. In the collision of term of the Boltzmann equation, the first term represents pair annihilations of ξ i andξ i into particles which are in thermal equilibrium (including the SM particles also) while the second term is due to the dark matter conversion process ξ 2ξ2 → ξ 1ξ1 where none of them are in thermal equilibrium during freeze-out. Detailed derivation of this term when both initial as well as final state particles are not in thermal contact with the visible world is given in appendix C. In the collision term, the thermal averaged annihilation cross section of a particular process AA → BB has been denoted by σv AA →BB . The possible annihilation channels of both the DM candidates in our model are shown in figure 2. This figure not only contains the Feynman diagrams for individual DM annihilations into other particles, but also the conversion of one particular DM pair into a pair of the other DM. We have solved these two coupled Boltzmann equations using micrOMEGAs [120] where the model information has been supplied to micrOMEGAs using FeynRules [121]. All the relevant annihilation cross sections of dark matter number changing processes required to solve the coupled equations are calculated using CalcHEP [122]. Finally, after solving the Boltzmann equations we get the comoving number density Y i (T 0 ) of each dark matter candidate at the present epoch (at T = T 0 ). Thereafter, one can easily calculate the total dark matter relic density which is sum of relic densities of all dark matter candidates, In order to understand how the comoving number densities are varying with respect to the temperature T , we have shown two plots in the both panels of figure 1 illustrating the thermal evolution of Y 1 and Y 2 for m ξ 2 > m ξ 1 (left panel) and m ξ 1 > m ξ 2 (right panel) respectively. In the left panel, we consider m ξ 2 = 2m ξ 1 with m ξ 1 = 2110 GeV. Here, the red solid line represents the variation of Y 1 with x ξ 1 = m ξ 1 T while the same for heavier component ξ 2 has been indicated by the blue solid line. The corresponding equilibrium number densities computed using the Maxwell-Boltzmann distribution are shown by red dash-dotted and blue dotted lines respectively. Form this plot it is seen that initially for low x ξ 1 (i.e. for high T ) the comoving number density of each component follows their respective equilibrium number density upto a certain temperature T (different for ξ 1 and ξ 2 ) and thereafter Y i departs significantly from the corresponding equilibrium distribution function Y eq i and remains constant with respect to the variation of T . This is nothing but the well known freeze-out point of a WIMP and it depends on when the interaction rate corresponding to the number changing processes of a particular dark matter component goes below the expansion rate of the Universe governed by the Hubble parameter H. In this particular situation, freeze-out of ξ 1 occurs at x ξ 1 ∼ 25 while that for the heavier component is . Furthermore, the interactions of ξ 1 and ξ 2 are such that ξ 2 which freezes-out earlier has less relic abundance compared to its lighter counter part ξ 1 . The opposite situation is presented in the right panel of figure 1, where ξ 1 is the heavier dark matter component. Here we consider the model parameters in such a way that although the heavier component ξ 1 has earlier freeze-out, ends up with a greater abundance.

Direct detection
Since each of the DM candidates in our model is a Dirac fermion, there exists Z BL as well as scalar mediated spin independent elastic scattering processes off nucleons. The relevant Feynman diagrams are shown in figure 3. Since several ongoing experiments like LUX [8], PandaX-II [9,10] and Xenon1T [11,12] are looking for such processes, regularly giving stringent upper bounds on DM-nucleon scattering cross section, we can further constrain our model parameters from these data. Parametrising DM and quark interactions with the dominant spin independent DM-nucleus scattering cross section can be written down as [123] σ SI where µ ξ i N is the reduced mass of DM nucleus system,b q are the quark-Z BL couplings which is same for all quarks in B − L gauge model. Also, A, Z are mass number and atomic number of the nucleus. From the Lagrangian of DM given in (3.14), we can write  L ⊃ [ξ i (λ ξ i s + λ ξ 1 p γ 5 )ξ i + q(λ qs + λ qp γ 5 )q]s j the corresponding spin-independent DM-nucleus scattering cross section can be written as Heref p,n are defined asf  gives f T G ≈ 0.91. We extract the spin independent elastic scattering cross section for both the DM candidates off nucleons from micrOMEGAs. Keeping the fact in mind that we are analysing a two-component DM scenario we have multiplied the elastic scattering crosssection by the relative number density of each DM candidate to find the individual effective DM-nucleon scattering cross section. As mentioned earlier, both gauge interactions as well as Yukawa interactions will contribute in the direct detection cross sections. However, the gauge mediated interaction is proportional to g 4 BL whereas the scalar mediated processes are mixing suppressed as well as Yukawa coupling suppressed (the Yukawa couplings of first generation of quarks are O(10 −6 )). Hence, in our present model DM scattering off nucleons through Z BL will have dominant contribution.

Results
As we have discussed in the previous section, this model predicts two automatically stable DM candidates ξ 1 and ξ 2 . The total relic density of DM can be expressed as, the sum of the relic densities of ξ 1 and ξ 2 , Ω DM h 2 = Ω ξ 1 h 2 + Ω ξ 2 h 2 where Ω ξ 1 h 2 and Ω ξ 2 h 2 are the relic abundances of ξ 1 and ξ 2 respectively. Let us now investigate the dependence of the relic densities, Ω ξ 1 and Ω ξ 2 , on the parameters of the model. In figure 4, we have shown the variation of relic density with the DM mass by assuming M ξ 2 = M ξ 1 . The other parameters of the model were chosen as M ψ 1 = 1.5 TeV, M ψ 2 = 2 TeV, M ψ 3 = 750 GeV, M s 1 = M s 2 = M s 3 = 1 TeV, M A 3 = 10 TeV, s 12 = s 13 = s 14 = s 24 = 0.2, M Z BL = 5 TeV, g BL = 0.3. The dashed (blue) and dotted (red) lines denote the relic densities of each DM candidate,Ω ξ 1 h 2 and Ω ξ 2 h 2 , whereas the solid (green) line is their sum (Ω DM h 2 ). The horizontal magenta line represents the relic density bound from PLANCK data [3]. Some important features of the model can be indicated here. The resonances due to the s-channel annihilation through Z BL and the different scalars (s 1 , s 2 , s 3 , A 2 , A 3 ) has reduced the relic density as expected.  annihilation into one scalar (M s i = 1 TeV) and one Z BL (M Z BL = 5 TeV) final state near 3 TeV mass. However, because of the resonance that effect is not visible. Another important point to note here is that, Ω ξ 2 is subdominant throughout the whole mass range. That can be explained as follows, ξ 2 has formed by combining N 2L and N 2R where as ξ 1 has formed from N 1L and N 1R . The B − L quantum number assigned for N 2L and N 2R is greater than the quantum number for N 1L and N 1R and that increases the annihilation cross section of ξ 2 by some numerical factor which results to smaller abundance.
One point to note here is that both the gauge coupling g BL as well as Yukawa couplings Y 1 , Y 2 play important roles in the DM annihilation. The gauge coupling g BL is involved in the annihilation channels mediated by gauge boson Z BL while the scalar mediated annihilation cross sections depend on Yukawa couplings Y 1 , Y 2 . In the lower mass region of DM g BL plays the main role whereas in the high mass region Y i also gives a significant contribution as Y i = √ 2M ξ i /u. In figure 5, we have shown the variation of the relic density for two different relations between M ξ 1 and M ξ 2 . In left panel, we have assumed M ξ 2 = 2 M ξ 1 and in right panel M ξ 2 = M ξ 1 /2 − other parameters remain same as in figure 4. In figure 5a, the resonances have occurred at different positions for Ω ξ 2 , at 250 GeV, 1.25 TeV, 1.63 TeV, and 2.5 TeV, whereas, Ω ξ 1 has the same behaviour as in 4 which is expected as we have assumed the M ξ 2 = 2 M ξ 1 . Figure 5b can be explained in a similar way.
In figure 6 we have shown the variation of the total relic density as a function of m ξ 1 for three benchmark values of g BL (0.01, 0.09, 0.4) and mixing angle (0.001, 0.01, 0.1). The left panel shows that the total DM abundance increases as we choose smaller values of gauge coupling. It is because small g will decrease the annihilation cross section and eventually increase the DM abundance. The right panel shows that the relic abundance hardly depends on the mixing angle. Here for simplicity we have assumed all the mixing angle to be same. In both the cases, we have assumed M ξ 2 = 2 M ξ 1 and the other parameters remain same as in figure 4.
As seen from the above plots, the relic abundance of both the DM candidates primarily depend upon the strength of their annihilation cross sections via scalar or gauge portal interactions. After understanding this behaviour from the plots with benchmark choices JHEP12(2019)109 of model parameters, we now move onto performing a random scan over a set of free parameters mentioned in table 3. While the physical masses of DM and all new particles in the model are varied in the mentioned range along with g BL , the singlet scalar mixing angles with SM Higgs boson as well as among themselves are kept fixed at 0.2. We then apply the relevant constraints one by one to arrive at the final allowed parameter space from all relevant bounds. We first show the parameter space in M ξ 1 − M ξ 2 plane allowed from total DM relic abundance, perturbativity, bounded from below criteria of the scalar potential in figure 7. In figure 8, we then show the parameter space in g BL − M Z BL plane allowed from total DM relic abundance, perturbativity, bounded from below criteria of the scalar potential. The bounds from the LEP and the LHC 6 are separately shown by the respective shaded regions which are excluded. The coloured bar in the two plots shown in figure 8 correspond to the relative abundances of two DM candidates ξ 1 (left panel) and ξ 2 (right panel) respectively.
In figure 9, we show the spin independent DM-nucleon scattering cross section for individual DM candidates as functions of their mass. All the points satisfy the total DM relic, perturbativity, bounded from below criteria of the scalar potential. We have also shown the projected sensitivities of XENONnT (blue region) and DARWIN (grey region) experiments in both the plots which clearly show that these two future experiments can probe a large region of parameter space of the model. As mentioned earlier, the actual scattering cross section is multiplied by individual relative number densities n ξ 1,2 /(n ξ 1 + n ξ 2 ) in order to compare with the XENON1T bounds [12] derived for single DM component. This clearly shows that the model remains very much sensitive to the direct detection experiments with many parts of parameter space already being ruled out. We then superimpose the collider  Table 3. The parameters of our model and ranges used in the random scan bounds also on these points and show the resulting parameter space in figure 10. The severe impact of XENON1T constraints is clearly visible here with a very few points allowed compared to the points allowed without applying the XENON1T bounds. Although the individual DM-nucleon scattering bound allowed several points in the parameter space to survive, as seen from figure 9, when we impose the condition that both the DM candidates must satisfy the XENON1T bounds, it results in a much smaller allowed parameter space, seen in figure 10. This may be due to the fact that the elastic scatterings of both dark matter components (ξ 1 and ξ 2 ) with the detector nuclei occur predominantly through Z BL exchange and hence, σ SI other hand, Z BL as well as scalar bosons mediated annihilation processes have significant impact in the relic densities of ξ 1 and ξ 2 . Therefore, in the latter case we have more parameters (Yukawa couplings, g BL etc.) and consequently the large portion of g BL − M Z BL plane is allowed. The allowed points in the chosen range of random scan will face further constraints from future LHC runs as well as direct detection experiments like LZ [13], XENONnT [14], DARWIN [15] and PandaX-30T [16].

LHC signatures of fermion triplets
Although typical LHC signatures of a U(1) B−L model is via search for dilepton resonance mediated by Z BL mentioned before, the present model can have additional prospects of being discovered at the LHC due to the presence of triplet fermions. While the neutral components of fermion triplets play the role of generating light neutrino masses through type III seesaw mechanism, the charged components can leave interesting signatures at colliders. Although a detailed study of collider aspects of such fermion triplets is beyond the scope of this present work and can be found elsewhere, we briefly highlight the additional advantage of having such triplets in a gauged B − L model. Detailed collider studies of charged components of fermion triplet can be found in the context of long lived wino in supersymmetric models; for example, see [124] and references therein. In the context of non-supersymmetric models, collider studies of such long lived fermion triplets can be found in [125,126] and references therein. We consider the third fermion triplet and its charged component as it can have enhanced production cross section at colliders due to larger B −L charge. It should be noted that unlike in usual type III seesaw model, here the charged components of fermion triplets can be produced at the LHC via both Z and Z BL gauge bosons, apart from usual photon mediation. The charged components can then decay into (i) charged lepton and Higgs boson or (ii) neutral fermion triplet plus on or off-shell W boson depending upon the mass splitting ∆M 80 GeV or ∆M 80 GeV. Although the components of a particular fermion triplet have degenerate masses (M k0 ) at tree level, one can make the charged components ψ ± k heavier by considering one loop electroweak radiative corrections [127,128], which result in a mass splitting ∆M ∼ 166 MeV between M ψ ± k and M ψ 0 k for M k0 1 TeV. Thus, in the second possible decay channel mentioned above, only the off-shell W boson is possible. The first decay mode is typically sub-dominant by the constraints from light neutrino masses, which require TeV scale fermion triplet Yukawa couplings with the SM leptons to be as small as around 10 −6 − 10 −5 . Therefore we consider the second decay mode only where the charged component of triplet decays into the neutral component and an off-shell W boson. Due to the small mass splitting, the dominant decay mode has final states ψ 0 k , π ± . The decay width of ψ ± 3 to ψ 0 3 and π ± can be written as 166 MeV, g is the SU(2) L gauge coupling, f π = 131 MeV [127] is the pion decay constant, and the first diagonal element of the CKM matrix V ud 0.974 respectively. Such tiny decay width keeps the lifetime of ψ ± 3 considerably long enough so that it can reach the detector before decaying. In fact, the ATLAS experiment at the LHC has already searched for such long-lived charged particles with lifetime ranging from 10 ps to 10 ns, with maximum sensitivity around 1 ns [129]. In the decay ψ ± 3 → ψ 0 3 π ± , the final state pion typically has very low momentum and it is not reconstructed in the detector.  On the other hand ψ 0 3 is a long lived particle as it can decay to SM leptons through the small mixing with the other triplet fermions and leaves the detector without interacting. Therefore, it gives rise to a signature where a charged particle leaves a track in the inner parts of the detector and then disappears leaving no tracks in the portions of the detector at higher radii.
Another crucial difference from usual type III seesaw extension of the SM is that here the production cross section of fermion triplets gets enhanced. The presence of Z BL significantly increases the production cross-section of the charged ψ ± 3 pairs at the collider. This is due to the fact that the decay of the BSM neutral gauge boson to the ψ ± 3 now happens on-shell in contrast to models where this decay takes place off-shell via SM Z boson and photon. Also, as there is no negative interference between the Z and Z BL mediated charged ψ ± 3 production channels, hence the addition of new channel always improves the production cross-section. In order to illustrate this improvement, we first show the variation of production cross-section σ pp→ψ + 3 ψ − 3 in the left panel of figure 11 for two different choices of centre of mass energy. As can be seen, the improvement in production cross section is more significant in 100 TeV centre of mass energies of proton proton collisions. Once produced, ψ ± 3 can give rise to disappearing charged track signatures mentioned above. The ATLAS experiment at the LHC put constraints on such disappearing charged track signatures for a long lived chargino decaying into a pion and wino dark matter, which is shown as the solid black line in right panel of figure 11. Although ψ 0 3 is not the DM candidate in our model, it effectively behaves like one at the LHC due to its long life. It can be seen that the existing LHC constraint can already rule out ψ ± 3 masses below 500 GeV from its searches for disappearing charged tracks, keeping the parameter space considered in this study within near future sensitivity. Also, the bounds will get slightly tighter in our model due to enhanced production cross section. We leave a detailed study of such type III seesaw signatures at the LHC in the framework of a B − L gauge model.
It should be noted that here we have adopted the narrow width approximation for Z BL which implies Γ Z BL /M Z BL 1. To justify this, we consider three benchmark scenarios namely, (i) Z BL can decay into both the DM candidates as well as the fermion triplet, (ii) Z BL can decay into the fermion triplet but not to the DM candidates, (iii) Z BL can decay into both the DM candidates but not to the fermion triplet. The corresponding widths, considering M Z BL = 5 TeV , g BL = 0.3 are shown below in tables 4 and 5 respectively. As can be seen from these results, even for the maximum possible decay width of Z BL we can have Γ Z BL /M Z BL < 5%. As discussed in several earlier works including [130,131] and references therein, narrow width approximation remains valid for such small decay widths.

JHEP12(2019)109 9 Conclusion
We have proposed a gauged U(1) B−L version of type III seesaw model which naturally predicts a two component Dirac fermion dark matter scenario due to the requirements of anomaly cancellation. Unlike type I seesaw scenario in U(1) B−L model where three right handed singlet neutrinos with B − L charge −1 each leads to cancellation of all anomalies without any need of additional chiral fermions, type III seesaw implementation leads to additional anomalies due to the non-trivial SU(2) L structure of triplet fermions. We show that a gauged U(1) B−L model with three fermion triplets required for type III seesaw can be anomaly free due to the presence of two neutral Dirac fermions having fractional B − L charges. Both of these fermions are naturally stable due to a remnant Z 2 × Z 2 symmetry. We study the DM phenomenology of the two component DM scenario in the model after incorporating all relevant theoretical and experimental bounds. Among the theoretical bounds, the bounded from below criteria of the scalar potential plays a crucial role restricting the parameter space. The relic abundance of the DM candidates are primarily dictated by their annihilations into SM particles mediated by Z BL as well as additional singlet scalars. We find the parameter space allowed from total relic abundance criteria and then apply the bounds from direct detection experiments on individual DM candidates. We find that all these constraints tightly constrain the parameter space which we scan through in our study, leaving a small region which can be probed at experiments operational at different frontiers. Finally we comment upon one interesting prospect of probing this model through production of charged components of fermion triplets at colliders. The presence of additional neutral gauge boson Z BL can significantly enhance the production cross section of charged triplets at LHC or 100 TeV future proton proton collider compared to usual type III seesaw extension of standard model. The charged components can then decay into the neutral ones and a pion with a long decay lifetime, leaving a disappearing charged track signature. Comparing with the ATLAS bounds on such signatures, we find that such triplets with masses upto 500 GeV or so can already be ruled out. While the charged components of fermion triplets can be probed via such disappearing charged track signatures, the neutral components, if long lived enough could be probed at proposed future experiments like MATHUSLA [132]. These triplets, however, do not affect the DM phenomenology much apart from opening another DM annihilation channel to triplet fermions whenever allowed kinematically.

JHEP12(2019)109
A Diagonalisation of mass matrix of real scalars Here, we have denoted c ij ≡ cos θ ij and s ij ≡ sin θ ij respectively. Now, using the elements of the diagonalising matrix, the physical scalar fields can be expressed as a linear combinations of unphysical states as

B Diagonalisation of mass matrix of pseudo scalars
The pseudo scalar mass matrix with respect to the basis state 1 √ 2 (A 1 A 2 A 3 ) T is given by

JHEP12(2019)109
As expected, the pseudo scalar mass matrix is also real symmetric. Therefore, the matrix M ps can be diagonalised by a 3 × 3 orthogonal matrix and the corresponding eigenvalues are given below, (B.1) C Collision term of the Boltzmann equation for dark matter conversion process ξ 2ξ2 → ξ 1ξ1 In this section, we have derived the Boltzmann equation for an interaction which converts one type of dark matter into another i.e. ξ 2ξ2 → ξ 1ξ1 . In this process, both initial and final states particles are dark matter candidates which are not in thermal equilibrium during and after their freeze-out. The situation is completely different when each type of dark matter annihilates into the SM particles in thermal equilibrium. The collision term of the Boltzmann equations depends strongly on whether the final state particles in a scattering process are in thermal contact with thermal bath of the Universe or not. This actually modifies the form of the collision term which reduces to the known form when outgoing particles are in thermal equilibrium. The contribution of the process ξ 2ξ2 → ξ 1ξ1 to the collision term of ξ 2 is given by = − dΠ ξ 2 dΠξ 2 dΠ ξ 1 dΠξ 1 (2π) 4 δ 4 (p ξ 2 + pξ 2 − p ξ 1 − pξ 1 ) M ξ 2ξ2 →ξ 1ξ1 2 × f ξ 2 (| p ξ 2 | , T ) fξ 2 ( pξ 2 , T ) − f ξ 1 (| p ξ 1 | , T ) fξ 1 ( pξ 1 , T ) , (C.1) where, p i is the four momentum of species i while the corresponding three momentum and energy are denoted by p i and E i respectively. The phase space factor dΠ i ≡ g i d 3 p i (2 π) 3 E i with g i is the internal degrees of freedom of species i having distribution function f i (| p i | , T ). Here we have neglected the factors related to Pauli blocking for fermions (dark matter candidates) as those are not significant enough for the present context. Furthermore, M ξ 2ξ2 →ξ 1ξ1 2 is the Feynman amplitude square of the concerned process ξ 2ξ2 → ξ 1ξ1 and it is spin-averaged over both initial and final states particles. The above collision term can also be as − dΠ ξ 2 dΠξ 2 dΠ ξ 1 dΠξ 1 (2π) 4 δ 4 (p ξ 2 + pξ 2 − p ξ 1 − pξ 1 ) M ξ 2ξ2 →ξ 1ξ1 2 (C.2)

. (C.4)
Now, the quantity within the curly brackets is equal to 4 σv E ξ 2 Eξ 2 with σ being the annihilation cross section of ξ 2ξ2 → ξ 1ξ1 . Substituting this in the above we get − 1 n eq ξ 2 n eq ξ 2 where, the quantity within curly brackets is the thermal averaged cross section σv ξ 2ξ2 →ξ 1ξ1 . Therefore, in terms of the thermal averaged cross section, the collision term of ξ 2 due to dark matter conversion process is given by − σv ξ 2ξ2 →ξ 1ξ1 n ξ 2 nξ 2 − n eq ξ 2 n eq ξ 1 n eq ξ 2 n eq ξ 1 n ξ 1 nξ 1 . (C. 6) In this work, we have assumed that there is no asymmetry between the number densities of dark matter candidate ξ i and its antiparticleξ i and therefore, n (eq) ξ i = n (eq) ξ i . Hence, the contribution of ξ 2ξ2 → ξ 1ξ1 to the collision term of ξ 2 is given by − σv ξ 2ξ2 →ξ 1ξ1 n 2 ξ 2 − n eq ξ 2 n eq ξ 1 2 n 2 ξ 1 , (C.7) while the collision term of ξ 1 also has an additional contribution which is exactly identical to the above except the overall −ve sign will be replaced by a +ve sign as the dark matter conversion process ξ 2ξ2 → ξ 1ξ1 increases the number densities of both ξ 1 andξ 1 .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.