Diraxiogenesis

The family of Dirac Seesaw models offers an intriguing alternative explanation for the smallness of neutrino masses without necessarily requiring microscopic lepton number violation, when compared to the more familiar class of Majorana Seesaws. A global $\text{U}(1)_\text{D}$ symmetry, that is explicitly broken by a higher dimensional scalar operator, ensures that the right handed neutrino does not couple directly to the Standard Model like Higgs and an exact gauged or residual lepton number symmetry prohibits all Majorana masses. We demonstrate that all three Dirac Seesaws possess a Pseudo-Nambu-Goldstone boson associated with the $\text{U}(1)_\text{D}$ symmetry, that we call the Diraxion, whose cosmological dynamics have so far been left unexplored. Furthermore we illustrate that a Dirac-Leptogenesis version of the recently proposed Lepto-Axiogenesis scenario can be realized in this class of models, leading to a unified origin of the observed baryon asymmetry and dark matter relic abundance. Explaining only the baryon asymmetry can lead to potentially observable amounts of right handed neutrino dark radiation with $\Delta N_\text{eff.}\lesssim 0.028$. On the other hand, if we only fix the dark matter abundance via the kinetic misalignment mechanism, this set-up could lead to detectable signatures in proposed cosmic neutrino background experiments via decays of eV-scale Diraxions to neutrinos. Here there is no domain wall problem, since topological defects decay to a subleading fraction of relic Diraxions. A key ingredient of all Axiogenesis scenarios is the dynamics of relatively light scalar called the Saxion, that in our case has a mass at the GeV-scale and which might reveal itself in heavy meson decays or collider searches. Our setup predicts isocurvature perturbations in baryons, dark matter and dark radiation sourced by fluctuations of the Saxion.


Introduction
Condensates of scalar fields, either elementary or composite, play an important role in our understanding of fermion and gauge boson mass generation. Apart from the celebrated Higgs mechanism these condensates also allow for the possibility of explaining either the observed dark matter relic abundance or the matter anti-matter asymmetry. Some of the earliest proposals for baryogenesis focused on the dynamics of oscillating scalar bosons [1].
With the advent of inflationary cosmology [2] it was realized that the vacuum expectation value (vev) of such a field can undergo a large excursion during the quasi de-Sitter phase of the early universe [3]. The requirement for this is the existence of a very flat or almost vanishing scalar potential, which is naturally realized in supersymmetric field theories such as the MSSM. About 300 of these flat directions, that only receive their masses from supersymmetry breaking and higher dimensional effective operators or radiative corrections, are known [4,5]. This prompted the development of the Affleck-Dine mechanism [6], where the large effective vev generates the baryon asymmetry from a (potentially Planck scale) suppressed baryon number violating scalar interaction. To transfer this asymmetry from the scalar sector to the baryons typically involves decays of the condensate to thermal bath particles, which can be automatically included if the flat direction is responsible for inflation and the subsequent era of reheating (see e.g. [7] for a recent example).
While the aforementioned scenario adheres to well known Sakharov-conditions [8] relying on CPT-conservation, there exist a class of scenarios in which CPT is spontaneously broken in the plasma of the early universe [9,10]: The Pseudo-Nambu-Goldstone boson (PNGB) θ of an abelian symmetry, such as e.g. baryon number, oscillates in a cosine-potential. The initial field value of its oscillations can be at most 2π times the decay constant f of the field. One typically requires a rather large f to ensure that the underlying symmetry was broken before inflation and to suppress isocurvature perturbations, see e.g. [11]. Since the baryon asymmetry in these models is proportional to the oscillation frequency of the PNGBθ , which is related to its mass, one typically also needs masses around [9,10] or far above the electroweak scale [11] making these proposals hard to test in laboratory experiments. Consult also reference [12] for a review and [13] for a systematic treatment of this scenario. The charge from the oscillating PNGB is transferred to the fermions ψ in the plasma via a (schematic) derivative coupling (∂ µ θ)ψγ µ ψ, which reduces toθ ψγ 0 ψ for a homogeneous and isotropic PNGB field. Following from the fact that this only time-dependent term breaks Lorentz symmetry spontaneously in the early universe plasma, one can deduce that CPT is also spontaneously violated, hence the name of this scenario.
A new approach pioneered by [14] under the name of "Axiogenesis" essentially combines the two aforementioned proposals with the phenomenologically attractive QCD axion [15][16][17][18][19][20][21]. Here one can have large decay constants f a and an observable PNGB whose mass is inversely proportional to f a , making the scenario potentially testable. A large field excursion of the radial mode of the complex scalar that houses the PNGB angle, known as the Saxion 1 , can be used to convert its oscillatory motion in the early universe into a coherent rotation withθ ̸ = 0. The rotation in the axion direction couples to the QCD anomaly and induces chiral asymmetries for the quarks via the QCD sphaleron, which gets reprocessed into a B+L asymmetry via the electroweak sphaleron process. On top of that one can use the axion field velocity for a new scenario of dark matter known as kinetic misalignment [22,23] (see also [24]), unlike the conventional misalignment which works under the assumption of negligible velocity [25][26][27]. However this comes with the drawback of overproducing dark matter if one fixes the baryon asymmetry to its observed value [14]. This conclusion can be avoided if the electroweak phase transition is modified [14], one introduces additional sphalerons from e.g. a gauged SU(2) R [28], one takes the chiral plasma instability into account [29] or the axion possesses very large couplings to the weak anomaly [30]. The basic scenario can also work for more generic axion-like particles [31], very heavy QCD axions [32] or with multiple scalars [33]. Additionally there can be rich gravitational wave signatures [34][35][36][37][38][39]. Another way to make Axiogenesis viable is to directly produce a B-L asymmetry instead of B+L via additional processes such as lepton number violating scatterings mediated by heavy Majorana neutrinos, known under the name of "Lepto-axiogenesis" [40][41][42], or R-parity violating (RPV) supersymmetry [43] (see reference [44] for a review on RPV).
In this work we follow the Lepto-axiogenesis route and try to explain the origin of the PNBG and the effective B-L violation from the same source. To this end we abandon the traditional QCD axion and focus on the PNGBs associated with neutrino mass generation.
While the corresponding particle for Majorana neutrinos with the fitting name Majoron [45] has been widely studied [46], we choose to focus on the case of Dirac neutrinos. Here there are typically two symmetries [47]: We assume global U(1) D to prevent the coupling between the Standard Model Higgs doublet and the right handed neutrinos together with a gauged or residual lepton number (or B-L) symmetry ensuring the absence of all Majornana masses. Non perturbative quantum gravitational effects manifesting as Planck-scale suppressed higher dimensional effective operators explicitly break the U(1) D symmetry and we call the associated PNGB Diraxion. The idea is that we produce equal and opposite asymmetries in the Standard Model fermions and the right handed singlet neutrinos due to the conservation of the total B-L, which never get equilibrated [48]. The electroweak sphaleron process is only sensitive to the the lepton doublet so effectively (B-L) SM is violated in the plasma and can be converted into a baryon asymmetry. A similar idea albeit with a QCD axion has been persued in [49] for the case of composite right handed neutrinos. However due to the unspecified UV nature of the non-perturbative mechanism behind the right handed neutrino formation the authors can only use the temperature from when on (B-L) SM is conserved as a free parameter. Our scenario relies on calculable, perturbative models realizing parametrically small Dirac neutrino masses from threshold corrections in the Seesaw spirit. We discuss the three possible cases for these Dirac Seesaws and show that the singlet scalar that appears in all those constructions can house the Saxion and Diraxion. Similar to the original spontaneous baryogenesis proposal [10], we will assume that the underlying global symmetry is only broken by a higher dimensional operator responsible for both the Diraxion mass and converting the Saxion oscillation into a Diraxion rotation. Our scenario is inspired by models combining scalar field Leptogenesis with neutrino mass generation such as spontaneous Leptogenesis from a very heavy oscillating Majoron [50] or inflationary Affleck-Dine Leptogenesis from the Type II Seesaw [51,52]. Affleck-Dine Dirac Leptogenesis from supersymmetric Dirac neutrinos involves two complex scalar fields, where the smallness of the F -term potential is connected to the tiny neutrino mass scale, and this scenario was investigated in [53,54]. Unlike these models we predict sub-eV PNGBs and are able to explain the dark matter relic abundance via kinetic misalignment [22] or parametric resonance [55,56]. If our setup is only responsible for baryogenesis, dark radiation in the ν R component can be produced with an amount of ∆N eff. ≲ 0.028. Alternatively if we reproduce (a fraction of) the dark matter relic abundance, the Diraxion can be heavy and metastable enough for decays to neutrinos, leading to a discovery potential in next generation experiments investigating the cosmic neutrino background. Requiring the cogenesis of the baryon asymmetry together with dark matter on the other hand leads to unobservably small ∆N eff. and too light or long-lived Diraxions. While isocurvature fluctuations in baryons and dark matter are a generic prediction of this kind of scenario, we find that our setup also produces dark radiation isocurvature modes and all three kinds of perturbations are induced by Saxion fluctutations during inflation.
Since the Diraxion only has suppressed two-loop couplings to photons we do not expect any signal in axion haloscopes or helioscopes. In our scenario we find slightly heavier Saxions than for Lepto-Axiogenesis [40] with masses around the GeV scale, which can be tested in collider experiments and rare decays of heavier mesons.
The relevant parameter space of our analysis is spanned by the Diraxion mass m a , the Saxion mass m S , the Diraxion decay constant today f a , the initial Saxion field value S i as well as the domain wall number N , which coincides with the mass dimension of the non-renormalizable operator responsible for the Diraxion mass. We can eliminate f a by requiring that the interaction of the Saxion with the thermal bath is slow, which is equivalent to fixing the amount of right handed neutrino dark radiation produced via Freeze-In. Successful Saxion thermalization from the Higgs portal interaction before an era of Saxion domination is viable for S i ≲ 0.1M Pl. . In order to keep the Diraxion light enough we fix N = 6. This allows us to predict the masses m a and m S for successful cogenesis of the baryon asymmetry and dark matter relic abundance.
Section 2 gives an overview over the family of Dirac Seesaws together with the details about the Saxion and Diraxion. We illustrate the cosmological evolution leading to Saxion oscillations and ultimately a coherent rotation in the Diraxion direction in section 3. Dirac-Lepto-Axiogenesis is the subject of section 4. Dark Matter and Dark Radiation are the topics of sections 5 and 6 respectively. Estimates of the Dark Matter and Dark Radiation isocurvature perturbations can be found in 7. Section 8 describes the required Saxion thermalization to avoid overclosure and excessive amounts of Dark Radiation. We discuss our findings for the regions of parameter space that realize both Leptogenesis and dark matter in 9 before we conclude and summarize in 10.

Dirac Weinberg-operator
The famous dimension five Weinberg-operator (LH)(L c iσ 2 H) withH ≡ iσ 2 H † is the only gauge invariant combination of Standard Model (SM) fields that generates a Majorana mass for the left chiral neutrinos [57]. However it has long been known that by adding gauge singlet right chiral neutrinos ν R together with a singlet scalar σ one can realize a model field SU ( analogue of the Weinberg-operator for Dirac neutrinos where c ν is a dimensionless Wilson-coefficient that appears together with a UV cut-offscale Λ UV ≫ v H and v H ≡ 246 GeV is the vacuum expectation value (vev) of the neutral component of the SM Higgs field. After σ condenses with the vev v σ the above operator leads to neutrino masses of To ensure the absence of any perturbatively (in field theory) or non-perturbatively (from quantum gravity) generated [58,59] Majorana mass terms for the SM neutrinos or heavy messenger fermions, we assume either an unbroken global symmetry like lepton number as a residual symmetry from a larger gauge symmetry, or a gauged symmetry like U(1) B-L , which can give rise to Dirac masses depending on the choice of charges and the breaking pattern [60]. In the following sections we discuss the tree-level UV-completions of the operator in (2.1). We will not concern ourselves too much with the details of the aforementioned symmetry for the absence of Majorana masses such as anomaly cancellation. Instead we focus on the cosmological evolution of the global U(1) D breaking. The explicit breaking Figure 1. Diagrammatic representation of the dimension 5 operators for the Type I Dirac Seesaw (left) and Type II Dirac Seesaw (right) giving rise to Dirac masses for the active neutrinos. For the Type II case one could also consider an additional insertion of σ not depicted here.
via quantum gravitational effects might induce the following higher dimensional operator contributing to the Dirac neutrino masses The largest contribution for v σ ≪ M Pl. comes from n = 1 and it is smaller than the contribution from the Dirac Weinberg-operator in (2.1) as long as

The three Dirac Seesaws
Here we introduce the three tree-level UV-completions to the Weinberg-operator in (2.1) and focus on the heavy new messenger fields whose threshold corrections lead to the tiny observed active neutrino masses. A systematic study of tree-and one-loop-level UV completions can be found in [61][62][63].

Type I
The most basic Dirac Seesaw [47,64] was discovered shortly after the more well known Majorana-Seesaw [65][66][67][68][69] and its phenomenology was explored in [70]. The mechanism requires nothing more than at least two generations of super-heavy, electrically neutral vector-like fermions N L,R and the field content and charges can be found in table 1.
We assume that their masses M N are not connected to U(1) D and a global or gauged U(1) B-L to be responsible for the absence of any Majorana masses. A diagrammatic representation of the mass generation mechanism was depicted on the left of figure 1. After integrating N L,R out, the light neutrino masses in the single flavor approximation read where we chose M N around the grand unification scale for illustration. Keeping M N fixed we can accommodate larger v σ by making the Yukawa couplings Y L Y R smaller. The oneloop correction to the active neutrino masses comes from closing the Higgs-singlet loop via an insertion of the mixed quartic coupling λ σH [47] δm and it is subleading for perturbative values of λ σH . This variant of the Seesaw can be embedded in the grand unified theories based on SU(5) [47] or SO(10) [71]. For example in the SO(10) case this comes with the drawback of spoiling matter unification: N R fills up the 16-dimensional spinorial representation together with the rest of the SM fermions and ν R , N L have to be added as additional gauge singlets. A simpler UV-completion based on U(1) B-L consists of giving all leptons L, ν R , N L , N R vector-like charges normalized to one implying that B-L can be gauged and then breaking it via the vev of a scalar φ with charge |Q B-L [φ]| > 2 [60], which has no direct couplings to fermions. This automatically forbids all renormalizable and effective operators leading to Majorana masses.

Type II
Instead of fermionic messengers fields that mix with the active neutrinos, for a Type II Seesaw one instead tilts a scalar potential slightly, to generate a tiny vev for a much heavier scalar field [72][73][74][75][76]. The first version of this mechanism applied to Dirac neutrinos was considered by [77] in the context of a two-Higgs-doublet model with a soft mass term µ 2 η † H, where we introduced a new doublet η that is a copy of the SM like doublet, but charged under U(1) D so that it couples to ν R This scenario was UV-completed in [78][79][80][81] by adding a singlet scalar σ whose vev generates the soft mass term: The scalar potential can contain the following two terms where we fixed Q D [σ] = 1 in both cases and one can observe that the scalar potential is identical to the two possible choices for the DFSZ axion model [20,21]. Note that these terms do not induce a mass for the imaginary component of σ, which becomes evident after diagonalizing the mass matrix for all pseudoscalars [79]. In the following we focus on the term ∝ κ which is linear in σ, because this allows us to discuss all Dirac Seesaw variants via the same effective operator (2.1). All required fields and charges can be found in table Figure 2. Diagrammatic representation of the dimension 5 operators for the Type III-a Dirac Seesaw (left) involving new triplet fermions and Type III-b Dirac Seesaw (right) involving new doublet fermions, both giving rise to Dirac masses for the active neutrinos.
1. If we assume that the positive µ 2 η is the largest scale in the scalar potential then the aforementioned trilinear term will induce a tiny vev that lies below the electroweak scale for κv σ ≪ µ 2 η and can naturally accommodate the observed neutrino mass scale without needing a small Yukawa coupling Y ν . A corresponding Feyman diagram can be found on the right hand side of 1. The Type II Seesaw scheme comes with an additional suppression factor κ/µ η when compared to the fermionic Type I Seesaw (for M N ≃ µ η ), which is why for the same v σ the neutral component of the η-doublet can be made lighter than the vectorlike neutrinos N . An estimate illustrates the interplay of the various scales Again the one-loop correction to the active neutrino masses comes from closing the Higgssinglet loop via an insertion of the mixed quartic coupling λ σH and it is subleading for perturbative values of λ σH . This scenario is the most straightforward to further UV-complete in terms of a gauged U(1) B-L since all gravitational and cubic anomalies will vanish, if one adds no other fermions than three generation of ν R with the same canonical lepton number Q B-L [L] = −1 as L. Renormalizable and effective Majorana masses for L, ν R can then be forbidden if the scalar φ dominantly responsible for the U(1) B-L -breaking has a charge |Q B-L [φ]| > 2 [60].

Type III
The Type III Dirac Seesaw scenario is defined by the presence of a hypercharge-less scalar iso-triplet ∆. In the first case, which we call Type III-a, ∆ has no couplings to SM fields and only interacts with ν R and the right chiral component of a vector-like fermion triplet T L,R [82] (see also [83]) (2.17) which was shown on the left of figure 2. Alternatively for the Type III-b model depicted on the right of 2, one may consider a direct coupling of ∆ to L by introducing vector-like doublet leptons D L,R [84] For both models all fields and charges where compiled in 1. Since an iso-triplet scalar will spoil the custodial symmetry of the SM scalar potential, its vev modifies the SM ρ-parameter [85,86] to be [82,84] The observed masses of the electroweak gauge bosons force ρ ≡ m 2 W /(m 2 Z cos(θ W ) 2 ) to be close to one, which requires v ∆ at or below the GeV-scale [87]. Indeed one of the motivations for this scenario is the observed deviation in the W -boson mass reported by the CDF collaboration [88], which can be stated in terms of the electroweak precision observable [89,90] known as the T -parameter [84] T The CDF-tension can then be explained by [82] v ∆ ≃ 4 GeV. (2.21) Such a low vev typically implies that the electrically neutral component of ∆ has to be below the weak scale as well, which can phenomenologically challenging 2 similarly to the light scalar in neutrino-philic Two-Higgs-doublet models [91][92][93][94][95]. On top of that the pesudoscalar component of ∆ playing the role of the Diraxion can typically have anomalous couplings to the weak SU(2) L gauge bosons, implying a different phenomenology compared to the case, where the Diraxion comes from an SM singlet. To avoid these problems we add an additional scalar singlet σ. An attractive way to generate a small v ∆ while keeping the components of ∆ ultra-heavy is the inclusion of a Type II suppression for v ∆ , completely analogous to the mechanism for a small v η in the previous paragraph 2.2.2. The relevant scalar potential reads and one finds [84] v which is suppressed compared to the electroweak scale as long as (2.24) Thus these versions of the Type III Dirac Seesaw correspond to the class of "nested Seesaws" [81,[96][97][98][99] with a built in double suppression from the heavy T or D masses and the tiny value . (2.25) The heavy fermion masses M F can be estimated from and as expected the vector-like iso-mulitplet fermions F = D, T are lighter than the singlet fermions N from the Type I Dirac Seesaw due to the additional suppression from the scalar sector. If we make no assumption about the value of v ∆ we find that the observed neutrino mass scale would require M F ≃ µ ∆ ≃ 10 8 GeV for order one couplings Y L , Y R , λ 4 and v σ = 10 6 GeV. The one loop correction is found from closing the Higgs-loop via the λ H coupling or by using the Higgs-singlet loop from the mixed quartic λ Hσ . We obtain a finite contribution of which is below the tree-level mass for perturbative values of λ Hσ . A straightforward way to UV complete these models in terms of a gauged U(1) B-L is to assume that the heavy fermions are vector-like under both the SM gauge group as well as B-L.

Saxion
The singlet scalar field has a renormalizable potential where we minimized the potential by balancing the tachyonic mass squared −|µ σ | 2 < 0 against the quartic coupling λ σ to define the nontrivial vev v σ ≡ |µ σ |/ √ 2λ σ . One can decompose the scalar field as where we anticipated that during inflation the radial mode will be displaced from its true vacuum Apart from the SM and ν R the particle spectrum will contain a massive scalar state S called the Saxion, whose mass squared reads The Saxion couplings were summarized in appendix B.

Diraxion
The most important ingredient of our scenario is the pseudoscalar Nambu-Goldstone-Boson (NGB) a of the U(1) D symmetry, also known as the Diron [47] or Diracon [79,80], which can be understood of the Dirac equivalent of the well-known Majoron [45,46]. Quantum effects like instantons or classical explicit breaking of the underlying global symmetry can generate a mass for a turning it into a Pseudo-NGB (PNGB), provided that its mass is parametrically below v σ and m S . In analogy to arguably the most well studied PNGB, the QCD-axion [15][16][17][18][19][20][21], we will call this CP-odd scalar the "Diraxion". This PNGB will not couple to any SM fermions and does not have any anomalies with the SM gauge group. Global symmetries like our U(1) D are expected to be broken by non-perturbative quantum gravitational effects [100][101][102], as can be seen from wormhole arguments. Heuristically these effects are encoded in Planck-scale suppressed non-renormalizable operators, that explicitly violate a given symmetry. Furthermore, if we also want to use the same operator to induce a "kick" in the angular Diraxion direction from the Saxion oscillation, we need an operator with dimension larger than five [40]. While we normalized the number of units by which the vev of σ breaks the global symmetry to be Q D [σ] = 1, the explicit breaking will in general occur with a different number of units and consequently we have to deal with a degenerate vacuum. The most straight-forward way to see this, is to consider a potential of the form which implies that the Diraxion decay constant is not v σ but instead 33) where N is the vacuum degeneracy factor also known as the Domain-Wall number and we define θ ≡ a f a . (2.34) One can deduce that the above interaction is invariant under an unbroken residual discrete Z N symmetry, where σ transforms as ω with ω N = 1. This mismatch between spontaneous and explicit breaking will have considerable implications for cosmology since it can lead to long-lived domain walls that might overclose the universe [103,104]. However domain walls can be attached to cosmic strings from e.g. the U(1) D breaking and the resulting hybrid defect will be unstable for N = 1 [105,106]. Furthermore even for N > 1 one can use the same effective operator responsible for the Diraxion mass to induce domain wall decays to Diraxions [107]. We begin with the simplest effective operator imaginable with more than five singlet fields Pl.
and one can see that the Diraxion mass will automatically be suppressed compared to f a without needing a small dimensionless coefficient c d as long as f a ≪ M Pl. and d > 4. This approach comes with the downside of having a potentially large number of domain walls N = d > 5. In section 5.3 it will be explained why our setup is also viable for domain wall numbers larger than one. Hence for simplicity we will only consider the operator in equation (2.35) for the rest of this paper as it relates two of the free parameters N and d. More operators and a way to generate a specific value of N can be found in appendix C. One can estimate e.g. for N = 6 that (2.36) however it turns out that we need a Diraxion mass between 1 meV and 1 eV, which would require a tiny coefficient c 6 . An exponential suppression of the Wilson-coefficient due to a large wormhole action [108,109] might accomplish this. Alternatively we illustrate in appendix C how the required effective operator could arise from a second scalar field φ, that we assume is charged under a gauged U(1) B-L , which is responsible for the absence of Majorana masses. One such operator is Pl.
+ h.c. (2.37) and the required charges read Here we assume that φ has no direct couplings to the Seesaw messenger fields. For the Diraxion mass we find in this scenario that (2.38) If we also charge σ under B-L then we have to use more complicated chiral charge assignments than just ±1 for the fermion fields in order forbid Majorana masses and and to avoid gauge anomalies. Since this will typically involve additional anomaly cancelling fermions, we refrain from writing down the explicit UV-completions.

Diraxion couplings
The Diraxion only has tree-level derivative couplings to ν R as well as the heavy BSM fermions and scalars. Here we focus on the low energy couplings to SM particles and ν R .

Derivative Couplings
From the Saxion kinetic term we find that the Diraxion-Saxion coupling reads We then remove the phase of σ by a field definition of the fields charged under U(1) D , which leads to derivative couplings from their kinetic terms. At low energies only ν R is in the plasma so we focus on its couplings which read (2.41) The couplings for the Type I and Type III-a Seesaw are reduced [31], because here the For the Type II Seesaw there are no additional fermions for ν R to mix with and the D R for a Type III-b Seesaw have the same charge Q D [T R ] = −1 as ν R . However for Type III-b there is also the mass mixing between D R and L with Q D [L] = 0 that induces a derivative coupling c ν L = α 2 for L. Such a coupling of the Diraxion to a fermion with hypercharge might be interesting as it can lead to chiral hypermagnetic instabilities [29]. It is evident that the parameter α ≪ 1 due to the heavy messenger fermions, which is why in practise we will take c ν R ≃ 1/N for all Seesaws and c ν L ≃ 0 for Type III-b. The effective coupling for the production of a single Diraxion from a neutrino line is given by Limits from SN1987A on the energy loss and deleptonization exclude 10 −12 ≲ g aν ≲ 10 −5 [110][111][112] and for larger couplings the Diraxion would remain trapped inside the supernova. Meson decays exclude only g aν ≳ 10 −3 [113]. It is obvious from (2.42) that our setup is compatible with these bounds.

Loop induced couplings
The Diraxion only couples to SM quarks and leptons via the one-loop diagrams in 3 and here we recast the results of [114,115] obtained for a Majoron in a conventional Type I Seesaw. For us the pseudo-scalar coupling to electrons is the most relevant and it turns out to be Stellar cooling arguments from the sun as well as red giants exclude g ae ≳ 10 −13 [116][117][118][119], which is not in conflict with our scenario. By closing the charged fermion line and attaching two photons to the one-loop diagrams depicted in 3, we obtain the two-loop Diraxion-photon coupling. Since our choice of U(1) D symmetry has no mixed anomalies with the SM gauge sector, the resulting coupling will be proportional to m 2 a /m 2 ψ [114], where m ψ is the mass of the fermion running in the loop. Consequently the dominant contribution comes from the electrons being the lightest SM fermion and it reads [114,115] (2. 44) In contrast to the case of a QCD axion the above has no contribution from mixing with the pions, as the Diraxion has no tree-level coupling to quarks. Due to their dependence on the tiny neutrino masses squared we find that these couplings are below all relevant constraints.
3 Two-field dynamics

Initial Saxion field value
We consider a radial mode whose vev is displaced during inflation from its minimum as S i ≫ N f a and treat S i as the intial condition for the evolution of the Saxion. Here we present two ways to induce the displacement of the Saxion field value during inflation, and in appendix D we also discuss a way to avoid or relax the isocurvature constraints. In the following we take the initial field value S i as a free parameter and m S (S i ) 2 ≃ 3λ σ S 2 i denotes the effective mass squared for a quartic potential where we defined the Saxion mass squared in our vaccum today to be m 2 S . The comoplex singlet scalar is assumed to be a spectator field during inflation.

Quantum fluctuations
It has been long known [120][121][122][123][124][125][126][127][128][129] that scalar fields with very shallow potentials undergo large field excursions during inflation as a consequence of quantum fluctuations. We can take this potentially large field value as an initial condition S i for the Saxion field. When the Hubble rate during inflation H I is much larger than the effective mass, quantum fluctuations grow the expectation value of S. The field excursion of S can be determined from the variance ⟨S 2 ⟩ of the Starobinsky-Yokoyama distribution [3] where m S is the bare mass defined in (2.31). From this it becomes apparent, why a large field excursion S i demands a very flat potential, meaning a small m S (or equivalently H I ≫ m S (S i )). The abvove estimate assumes that only a single field is excited in this way [130]. It requires N e ≳ (H I /m S (S i )) 2 many e-folds of inflation for S to get pushed to the value S i [3]. One finds that the correlation of the associated fluctuations scales as [127] l ≃ H −1 and it is much larger than the horizon size H −1 I as long as H I ≫ m S (S i ), justifying the treatment in terms of a homogenous Saxion condensate.

Non-minimal coupling to gravity
While quantum fluctuations required m S (S i ) ≪ H I a non minimal coupling to gravity can lead to a Hubble dependent mass m S (S i ) ≃ H I . This kind of mass term was first suggested in the context of supersymmetry breaking from finite energy density effects [4,131] in the early universe. The same net effect can be generated by coupling the singlet scalar to the Ricci scalar . (3.5) The idea is to balance this tachyonic mass squared during inflation against the quartic term λ σ S 4 /4, which requires c R > 0. In order to have m S (S i ) 2 = −H 2 I we take c R = 1/12 and find The Saxion remains stuck in this value even at the onset of radiation domination because of the equilibrium between the Hubble friction and the slope of the potential [133][134][135]. Before we close let us note that an epoch of intermediate matter domination after inflation with R = −3H 2 , such as e.g. the well motivated curvaton scenario [136][137][138], could also be used to induce the large field value, which would be entirely free from inflationary isocurvature perturbations. The absence of fluctuations during inflation requires and effective mass during inflation that is larger than H I , which could be realized via the couplings to the inflaton introduced in in appendix D.

Saxion oscillation
The Hubble rate during radiation domination (RD)t and non-instantaneous reheating (RH), which corresponds to matter domination, in terms of the number of relativistic degrees of freedom g * . The radial mode starts oscillating with a frequency m S (S i ) when 3H(T osc. ) = m S (S i ) which determines for radiation domination The oscillation occurs after reheating (T osc. < T RH ) as long as For a quartic potential during radiation domiantion the field value S ∼ 1/a redshifts like radiation as a function of temperature for S > N f a After the Saxion starts to relax to its minimum S = N f a at the temperature the potential is dominated by the constant mass term in (2.31) and the field value of the now non-relativistic condensate redshifts as (3.13) For oscillations before the completion of reheating we find T RH (3.14) instead. During reheating we find the scaling law S ∼ 1/a ∼ T 8/3 because during matter domination T ∼ 1/a 3/8 . For most of our parameter space T S is below T RH , so the estimate in (3.12) applies with T RD osc. replaced by T RH and S(T RH ) obtained from the previous scaling. We can also compute the number density of the oscillating and thereby non-relativistic Saxion condensate which follows from its energy density ρ S = V σ (S) = m S n S . We want the evolution of the cosmological background to proceed as usual, which is why we demand that there is no period of inflation due the large Saxion field value [139] V By making use of m S (S i ) = 3H(T osc. ) we find that this implies [23,37] The evolution of the Saxion field during radiation domination, its energy density and its equation of state parameter ω S can be summarized by the following scaling relations In this work we assume that the masses of all additional particles are given by their treelevel expressions and do not get modified due to the large initial field value S i . Specifically this means that If this condition was violated ν R and N L would form a Dirac fermion of mass Y R S i or H, η would mix strongly. Furthermore we impose that the heavy messenger fields are not present in the plasma at any point in time where H I /(2π) is the Gibbons-Hawking temperature during inflation [140], T RH the radiation bath temperature at the end of reheating and the maximum temperature during (non-instantaneous) reheating [141,142], which can be much larger than T RH . CMB observation constrain the Hubble rate during inflation to be H I < 6 × 10 13 GeV [143] and thus puts an upper limit on the energy scale during inflation from which one can deduce T RH < T max < 6.5 × 10 15 GeV [29]. Note that it will in general be hard to satisfy (3.21) for the Type III scenario due to the smallness of µ ∆ (see equation (2.23)) unless we assume a low Hubble scale during inflation and a small reheating temperature. As long as (3.21) is satisfied the Saxion will not receive a thermal mass of the order Y T , where Y is a Yukawa or the square root of quartic coupling from (3.20). However since the heavy η or ∆ coupling to σ have electroweak gauge interactions, integrating them out can potentially modify the logarithmic running of the weak gauge couplings manifesting itself in a so called thermal logarithmic potential [144,145] Here c W is a model-dependent coefficient from the couplings of η, ∆ to σ. Demanding that this correction is subdominant to the quartic potential, so that the Saxion oscillates around its mass in (3.1), leads to the condition .
To be conservative we will take S i > 10 16 GeV throughout this work. Note that this bound will be absent for the Type I Seesaw, because the heavy vectorlike neutrinos have no weak gauge couplings. In section 7.3 we investigate the one loop quantum corrections to the Saxion quartic.

Generating the Diraxion rotation
We assume that the global symmetry is spontaneously broken with a vev S i during inflation that is significantly larger than today. Due to this vev, Planck scale suppressed explicit symmetry breaking will be active. As longs as the mass of the angular mode (during inflation) is small compared to the radial mode, we can consider U(1) D as an approximate symmetry (during inflation). The explicit breaking then becomes irrelevant as the vev relaxes to its true vacuum at N f a today. During inflation the Diraxion will have a mass In order for our description in terms of an approximate U(1) D symmetry to apply we require that the Diraxion mass during inflation is smaller than the Saxion mass in (3.1) [32], which leads to For N = 4 this would reduce to the requirement for the PNGB mass compared to the mass for the Higgs scalar of the underlying symmetry m a < m S . We further deduce that we can not take N to be arbitrarily large in order not to spoil the previous relation. Next one defines the charge yield or asymmetry of the scalar condensate as which is conserved as long as there are no interactions that explicitly violate the underlying global U(1) D symmetry. One can understand the fact that only a rotating condensate can carry a charge compared to an oscillating one from the observation that an oscillation is a superposition of two rotations in opposite directions [146]. The Planck suppressed effective operators in V / D of section 2.4, that are responsible for generating the Diraxion mass, convert a part of the oscillatory motion of the radial mode into a velocity for the angle. After the radial mode has decreased via redshifting and eventually settled to its minimum, the Planck suppressed operators σ N become negligible for N > 5 [40].The corresponding equation of motion reads [40] Pl.
Here θ i is the initial angle of the Diraxion also known as the misalignment angle. It is customary (see e.g. [147,148]) to define the eccentricity parameter ε [14,40] where we made use of the definition of the Diraxion mass in (2.35) and the Saxion mass in (3.1). The last approximate equality holds for sin(θ i + δ) ≃ O(1), which suppresses axion isocurvature perturbations [40] (see section 7). This parameter ε ≤ 1 has the following physical interpretation: An orbit with ε = 1 corresponds to a perfectly circular rotation, whereas one would have ε = 0 for a pure oscillation in the Saxion direction. Our scenario is different from the case of QCD-Axiogenesis [14] because here the "kick" in the angular direction and the Diraxion mass originate from the same operator, so that we can in principle choose ε ≃ 1. For a quartic potential and a QCD axion this is in general not possible [23], because of the axion quality problem [149][150][151] from the Planck suppressed operators, which can shift the minimum of the QCD axion to too large angles compared to the limit from the neutron's electric dipole moment. However for ε = 1 two complications arise: On the one hand we can no longer neglect V / D compared to the quartic term when determining the Saxion mass and initial field value, so our present description based on (3.1) breaks down. Additionally for larger V / D one can no longer neglect the angular gradient compared to the Hubble friction, so that the Diraxion stops being overdamped and starts to relax to the minimum of its potential given by (see also the discussion in section 7.1) If the Diraxion is trapped in such a vacuum, the rotation of the condensate will not occur or be significantly damped [40], which is why we demand ε < 1 This condition intuitively means that the Diraxion should not start oscillating before the Saxion [40], which is nothing more than the requirement m a (S i ) < m S (S i ) we already imposed in (3.26) to have a Pseudo-Nambu-Goldstone boson, hence the bounds being the same up to a prefactor. Both conditions automatically make sure that Diraxion is in motion [32] since its velocityθ ≃ m S (S) is always larger than its mass m a (S). In a quartic potential the angular and radial modes redshift as radiation for S ≫ N f a (ω S = ω θ = 1/3). For cosmological reasons that will be elaborated upon in section 8, the Saxion needs to be thermalized and will typically loose its energy to the SM plasma. After the damping of the Saxion oscillations only the angular rotation, which is now perfectly circular, remains. The angular rotation is stable because of U(1) D -charge conservation. It minimizes the free energy [14] and can loose energy to thermal bath via interactions with particles charged under U(1) D leading to the production of fermionic asymmetries. However as long as the Noether charge of the condensateθS 2 is larger than the typically produced fermionic asymmetryθT 2 , it is energetically favoured for the rotating condensate to retain most of its charge compared to the production of fermions. This requirement can be expressed as [14,152] where T S defined in (3.12) is the temperature at which the Saxion reaches its minimum S = N f a . This condition implies that T osc. ≪ S i for a quartic potential with the oscillation temperature from (3.8) and it can be expressed as which is not a strong constraint compared to the limits from isocurvature in section 7.1. The tilt in the angular direction responsible for the Diraxion mass can also lead to washout by converting Diraxions into Saxions, who scatter of the plasma with a rate Γ S defined in sections 4.1 and 8. The resulting rate for the angular mode is then suppressed by a factor of m a which is too suppressed to matter compared to Γ S . Once it thermalizes and reaches S = N f a the Saxion behaves as non-relativistic matter (ω S = 0) and the Diraxion has the equation of state of kination ω θ = 1, because only the kinetic energy of the remaining circular rotation is left, which is why then ε = 1. We can summarize the evolution of the angular field by the following scaling relationṡ We further define the conserved charge yield n θ /s for an oscillation starting during radiation If the oscillation starts before the end of reheating T osc. > T RH we need a different adiabatic invariant than n θ /s, because entropy is not conserved during reheating. Instead one normalizes n θ ∼ 1/a 3 to the energy density of the non-relativistic inflaton (3.39) After Saxion thermalization described in section 8, where ε → 1, and while the Saxion still approaches its minimum N f a one finds that the equations of motion for θ reduce to the balance between the radial potential gradient and the centripetal force S(θ/N ) 2 [29] Note that the energy density in the rotation ρ θ = ερ S (with ρ S = V σ (S i )) does not affect the condition of having less energy in the condensate than the radiation bath leading to (3.17), because the Diraxion energy originates from the Saxion oscillation, whose energy decreases to (1 − ε)ρ S after the kick in the angular direction. The 1/a 6 scaling of the Diraxion's energy density after Saxion thermalization is reminiscent of kination. For a quartic potential starting out during radiation domination it was argued in [37] that no era of kination dominance can occur without prior intermediate matter domination: This follows from the fact that the thermalization of the radial mode produces a radiation bath of energy density (1 − ε)ρ S , which is comparable in size to the energy density in Diraxions ερ S . Even if the Diraxion contribution is initially larger than the radiation bath, it will redshift faster and quickly become subdominant.
As a consequence of the oscillations for ε ≪ 1 the quantities S andθ will not be constant in time and oscillate themselves around their minimum and maximum values during one cycle which can be obtained from the conservation of charge and energy where S max denotes the maximum field value during a cycle and we identify The Diraxion velocity is as large asθ max for a time scale ∆t ≃ 1/θ max . If we compute the cycle average ofθ over a time scale t S ≡ 1/m S (S) we find that [40] ⟨θ⟩ ≃ ∆t More refined analytical [37] and numerical [40] calculations have confirmed that ⟨θ⟩ = N m S (S) is indeed an attractor solution, meaning that ⟨θ⟩ is independent of ε in practise [40]. In the above expression for m S we slightly abuse notation by denoting the root-meansquare ⟨S 2 ⟩ of the Saxion amplitude with the same symbol as its (oscillating) field value S. If parametric resonance from Saxion oscillations occurs (see section 5.5) the average oḟ θ gets reduced and one finds [40] ⟨θ⟩ = εN m S (S). (3.45) Note that reference [40] employs a slightly different parameterization of ε This parameterization reduces to our choice of ε in (3.31) for ε ≪ 1. We assume thatθ and consequentlyθ max changes only adiabatically, meaning slower than the time scale of the plasma 1/T , which is why we imposeθ max < T osc. [40] leading to m a > 0.24 Additionally we assume that the Saxion field value is always the largest energy scale in the plasma S min = εS i > T osc. [40] from which we find and one observes that both bounds are comparable in strength. These bounds are weakest for N = 5 due to an approximate cancellation in the exponent of N f a /S i , and get stronger for larger N . For reference we will take N = 6 throughout this paper.

Dissipation coefficient from Dirac Weinberg operator
The Saxion only couples to the thermal bath via its mixed quartic with ther Higgs or the non-renormalizable Dirac Weinberg operator. A proper determination of the relevant dissipation coefficient [153,154] for our scenario would involve two-loop thermal self energy diagrams [155], which is why we resort to dimensional analysis. The Saxion can scatter via the reactions HS → Lν R , LS → H † ν R or decay as S → LH † ν R , where we assumed a vanishing thermal abundance of ν R and put the non-thermal S in the initial state. Due to time-dilation effects [156] the decay will only be relevant for m S (S) < T , which is why we add both contributions as The scaling is similar to the interaction of two Higgsinos with two Saxions considered in [42] and can be understood as follows: The matrix element involves one insertion of S and the factor m ν /(v H N f a ) from the coefficient of the Dirac Weinberg operator. After squaring the matrix element to obtain the rate one can use dimensional analysis to find the remaining factor of 2T for the two scattering modes or m S (S)) for the decay. One should keep in mind that the rate depends on the time-dependent field value S and not on its amplitude S max [157]. After thermalization of S the rate reduces to which is now slower since T ≪ S. The overall rate for T > m S (S) differs only by a factor of v 2 H /(N f a ) 2 from the the result for Majorana neutrinos. We always find that m S (S) ≪ T for S ≫ N f a because m S (S i ) = 3H(T osc. ) ≪ T osc. , which is why the Saxion decay is negligible above T S . After the Saxion has reached its minimum at T S defined in (3.12) and we replace it by its vev, we recover the Yukawa interaction between the neutrinos and the Higgs, which never thermalizes due to the smallness of m ν /v H [48].

Chemical potentials and Boltzmann equation
The derivative coupling of the spatially homogeneous Diraxion ∂ µ θν R γ µ ν R defined in (2.41) can be thought of as inducing an effective chemical potential for the right chiral neutrinos [9,10,158] As pointed out in [159][160][161] and more recently in [162] however this derivative coupling does not lead to actual splittings in the single particle energies of particles and anti-particles, which is why it is not an actual chemical potential. A more accurate way to think about the effect ofθ is to treat it as a background field that shifts the dispersion relation of the particles coupling to it E = p 2 + m 2 → p 2 + m 2 ∓ c ν Rθ [159][160][161] with a (minus) plus sign for (anti-)particles in analogy to explicit CPT-violation in the form of different masses for particles and antiparticles, see e.g. [163]. Since the background fieldθ ̸ = 0 spontaneously violates CPT in the plasma, this scenario is known as spontaneous baryogenesis [9, 10] and we do not need to worry about the third Sakharov condition [164], which assumes CPT-conservation. The time dependent Diraxion background is also responsible for the CP violation encoded in the second Sakharov condition [164]. If the Dirac Weinberg operator leads to reactions HS → Lν R , HL → Sν R or S → Lν R H † in equilibrium, then the chemical potentials satisfy where we used that the real valued Saxion has no chemical potential. The same result can be obtained without making reference to the effective chemical potential µ eff. ν R : Suppose we work in the basis, where we do not remove the phase of the Dirac Weinberg operator m ν /(v H N f a )e iθ(t) and hence have no derivative couplings of the Diraxion to ν R . The time dependent phase in the aforementioned coupling does not enter the matrix element but rather it modifies the delta distributions encoding energy-momentum-conservation by also taking the effect of the slowly varying background field into account [50]. For e.g. [159][160][161] and via the collision term of the Boltzmann equation 3 one is then lead to the same equation for the chemical potentials as in (4.4). We now also impose the conservation of the total B-L, which might be gauged in the UV and is responsible for the Dirac nature of neutrinos, Here the actual chemical potential µ ν R appears, because as previously discussed c ν Rθ only contributes to scattering processes involving ν R . Next, if we impose the conservation of hypercharge, take all SM Yukawa and gauge interactions to be fast (see e.g. appendix C in [165]) and treat all flavors the same, then the network of chemical reaction implies that in equilibrium The relation between µ eq. B and µ SM eq. B-L is given by the standard value 28/79 [166,167] for the sphaleron redistribution coefficient. Note that the Yukawa interactions for all three generations of SM fermions are only in equilibrium below about 10 8 GeV [168], but we do not expect this detail to change the sphaleron redistribution coefficient by more than 10%. By employing the relation n i = µ i /6T 2 between asymmetries and chemical potential for a fermionic particle i and using the principle of detailed balance we can immediately write down the Boltzmann equation for n SM in terms of the dissipation coefficient Γ S (S) defined in (4.1). Our mechanism is similar to the case of Dirac-Leptogenesis [48] in the sense that total lepton number vanishes and we produce equal asymmetries (n SM B-L = n ν R ) in the left-and right-chiral leptons. These leptonic asymmetries are never equilibrated because in the following we take Γ S to be slow for S ≫ N f a . Once S ≃ N f a we recover the usual Yukawa interaction with the Higgs from the Dirac Weinberg-operator, which is always slow due to the smallness of the effective Yukawa coupling m ν /v H [48]. The B+L-violating sphaleron transition only acts on the doublet of left-chiral leptons and converts n SM B-L into the observed baryon asymmetry. A sketch of the scenario can be found in figure 4. Conventional Lepto-and Baryogenesis from scattering was covered in [169][170][171].

Baryon Asymmetry
We rewrite (4.7) in terms of the dimensionless yield n SM B-L /s, which is conserved in the comoving volume as long as entropy is conserved Solving (4.8) is in general complicated by the fact that our problem has two intrinsic timescales: The expansion rate of the universe and the oscillation frequency m S (S). Since our dissipation coefficient is given in terms of an effective operator, it will be UV-dominated and hence we expect the baryon asymmetry to be predominantly produced at T osc. . By definition we have 3H(T osc. ) = m S (S i ) at this point in time, which simplifies the problem. The last complication is the fact that the production rate Γ S (S) depends on the oscillating field value, which is why we need to average over one cycle of oscillations. Reference [40] solved (4.8) numerically (see appendix E of the aforementioned reference) and provided analytical arguments to understand the resulting abundances given by . (4.9) The first result can be understood as follows [40]: In order for ⟨n SM B-L ⟩ /s to track its equilibrium value,θ needs to be close to its maximumθ max = m S (S i )/ε > m S (S i ). However due to the conservation of the condensate charge n θ = εm S (S i )S 2 i we find that this only occurs when S is at its minimum for a cycle given by S min = εS i (see (3.41)). For a fielddependent interaction it was found in [40], that the charge transfer from the condensate to the plasma is only efficient if the transfer rate is larger than the frequency of the coherent motion, which for this case reads and reduces to the condition Γ S (S i ) ≫ m S (S i )/ε 3 displayed in (4.9). However if we take e.g. ε = 0.1 then this implies that Γ S (S i ) would have to be a thousand times faster than the Hubble rate at the beginning of the oscillations. We discuss in section 8.1 how this is potentially problematic. The next regime is now too small to trackθ max , the transfer can only be relevant when S is close to its maximum value S max (T osc. ) = S i . Once S decreases from its maximum during a cycle andθ increases, ⟨n SM B-L /s⟩ stops trackingθT 2 /s as soon asθ reaches . (4.11)

Rotation during Radiation domination
To be conservative we focus on the regime Γ S (S i ) ≪ m S (S i ) ≃ 3H(T osc. ), which can be thought of as Freeze-In production [172] in contrast to the previous two cases, which are akin to thermal Freeze-Out. We average the Boltzmann equation in (4.8) over a time scale larger than t S = 1/m S (S i ) and factor out the S-dependence of the rate by employing (4.1) and (4.2) to rewrite Γ S = Γ L · (S/T ) 2 leading to d dt Because θ and S move together coherently in the broken phase, we have to average them together as ⟨θS 2 ⟩ (unless parametric resonance occurs).
and obtain ⟨n SM B-L /s⟩ ∼ ⟨θS 2 ⟩ / ⟨S 2 ⟩ = εN m S (S i ), which explains the factor of ε in the last line of (4.9) and the factor of two comes from matching to the previous two solutions.
However d dt n SM B-L s = 0 can also be obtained from Γ S → 0 and we expect [172] the Freeze-In abundance to vanish if the production rate goes to zero (there is only an upper limit on Γ S in the third line of (4.9)). Therefore we write where in the last line we used the charged yield during radiation domination defined in (3.38). The same result can also be found from dropping the back-reaction term ∝ ⟨n SM B-L /s⟩ in (4.12) and considering the amount of asymmetry produced ⟨n B-L /s⟩ ≃ 1/H · d/dt ⟨n B-L /s⟩ over one Hubble time 1/H together with the charge conservation in (4.13). One can observe that the resulting asymmetry is suppressed by both ε < 1 and Γ S /H ≪ 1 compared to the scenario for a very efficient charge transfer Γ S (S i ) ≫ m S (S i )/ε 3 , which also justifies neglecting the back-reaction. By using the scaling relations for the redshift of S one can deduce that during radiation domination we have m S (S) ∼ S ∼ T and Γ S /H ∼ S 2 /T ∼ T leading to ⟨n B /s⟩ ∼ T , which means that the produced asymmetry is indeed UV dominated. We deviate from [40] by writing our result in terms of m S (S i ) instead of ⟨θ⟩ for the following reason: Typically one has ⟨θ⟩ = m S (S i ) at the beginning of oscillations, except when parametric resonance from Saxion oscillations takes place (see section 5.5), where one obtains ⟨θ⟩ = εm S (S i ) instead. However our asymmetry is sourced by ⟨θS 2 ⟩ and not just ⟨θ⟩, so no additional factor of ε arises. The net effect of parametric resonance is that the radial and angular fields will not move coherently anymore due to the non-thermal symmetry restoration, which simply means that we can separate the averages ⟨θS 2 ⟩ = ⟨θ⟩ ⟨S 2 ⟩ [40]. 4 The baryon asymmetry differs from the result for Majorana Lepto-Axiogenesis [40] in two aspects: On the one hand our asymmetry always depends on ε due to the linear coupling to S and on the other hand the ratio Γ S /H is parametrically different from the suppression factor Γ L /H found in the Majorana scenario. The order of magnitude and sign of the asymmetry will be determined by the first oscillation [11,173,174]. If we plug in the out-of-equilibrium condition for successful Freeze-In Γ S (S i ) ≪ m S (S i ) = 3H(T osc. ) we obtain an upper limit on the produced asymmetry 16) in terms of the charge yield of the condensate. Since the above is much less 5 than n θ /s, it was valid to neglect the back-reaction of the asymmetry production on the condensate; Leptogenesis will not evaporate the condensate. The produced leptonic asymmetry is not washed out, because the Dirac Weinberg-operator is the only coupling connecting ν R to the bath and we take it to be slow throughout the evolution of the universe. Once the Saxion settles at its minimum N f a the Dirac Weinberg operator reduces to the Yukawa coupling of the neutrinos with the SM like Higgs, which never equilibrates [48]. To estimate the baryon asymmetry we eliminate f a via the condition Γ S (S i )/H(T osc ) < 1 The baryon asymmetry only depends on the ratio S i /f a , so the dependence on S i divides out when inserting (4.17). This parameter choice is compatible with the bound from ε < 1 in (3.33).

Rotation during Reheating
In the above we assumed that the oscillation occurs during radiation domination. For the opposite case of T osc. > T RH we have production during reheating, where entropy is not conserved. Thus n B /s is not the correct adiabatic invariant anymore and we have to use ⟨n B ⟩ /ρ inf evaluated at production multiplied with ρ inf. /s evaluated at reheating. In this context ρ inf denotes the energy density of the non-relativistic inflaton. During the matter dominated reheating phase we have T ∼ 1/a 3/8 and a Hubble rate of H ∼ T 4 . By using the Friedmann equation ρ inf. = 3M 2 Pl. H 2 we find that ρ inf. ∼ T 8 . The scaling of the asymmetry 4 Using this together with the charge conservation (4.13) then readily implies ⟨θ⟩ = εmS(Si). 5 This assumes that Si > Tosc/ √ N , which is valid throughout our parameter space and we used cν R ≃ 1/N . could be different depending on whether the Saxion settles to its minimum N f a during (T osc. > T S > T RH ) or after reheating (T osc. > T RH > T S ). For T RH > T S we find that m S (S) ∼ S ∼ 1/a ∼ T 8/3 and Γ S (S)/H ∼ S 2 /T 3 ∼ T 7/3 implying ⟨n B ⟩ /ρ inf. ∼ 1/T . In the second regime T RH < T S we have S = N f a = const. so conservation of the Noether chargeθN 2 f 2 a implies ⟨θ⟩ ∼ 1/a 3 ∼ T 8 . Further we find Γ S /H ∼ T f 2 a /T 4 ∼ 1/T 3 resulting in ⟨n B ⟩ /ρ inf. ∼ 1/T . Both regimes redshift the same, because they depend on ⟨θS 2 ⟩ ∼ n θ , which always redshifts as 1/a 3 no matter if S is at or above N f a . During reheating the asymmetry is always IR-dominated and gets predominantly generated at T RH , which is why we can evaluate ⟨n B /s⟩ at reheating. This is unlike the usual case for Lepto-Axiogensis, where one would find ⟨n B ⟩ /ρ inf. ∼ T for T S > T RH [40]. We arrive at where in the second line we used the charge yield for reheating defined in (3.39) as well as the scaling laws to relate the field value at reheating S RH and S i , which restores the factor of Γ S (S i )/H(T osc. ) ≪ 1 6 . Note that the above is not enhanced by · T RH 10 13 GeV

Dark Matter decay
The Diraxion is only a good dark matter candidate if it is sufficiently long-lived. Its main decay mode is to neutrinos, followed by a subleading mode to photons via a two-loop coupling and the width is given by [175] Γ (a → νν) ≃ where we estimated the phase space suppression in the single flavor approximation. Since the coupling is dependent on the absolute neutrino masses ν m 2 ν we take the lightest neutrino mass to be zero so that where NH (IH) denotes the normal (inverted) hierarchy of neutrino masses. Throughout this work we focus on the normal hierarchy. Normalized to the age of the universe t 0 ≃ 13.79 Gyr we obtain Experiments such as PTOLEMY [176,177] aimed at detecting the cosmic neutrino background could potentially detect [178,179] the neutrinos produced in decays of dark matter with masses between 0.1 eV and O(1 eV) for lifetimes between (10 − 100) t 0 [175]. We show the corresponding parameter region colored in red in figure 5. While standard misalignment and decays of topological defects (blue region in 5) works only for too large values of f a compared to the detectable region, we will see that kinetic misalignment (see the orange lines in the aforementioned plot together with 5.4) and parametric resonance in 5.5 could both reproduce the relic abundance for much smaller f a and could potentially be probed by PTOLEMY.

Standard Misalignment
First we discuss the conventional misalignment mechanism, which operates when the kinetic energy of the Diraxion is smaller than its potential barrier. We will state the precise conditions for this at the end of this subsection and in the next one. During inflation the Diraxion is trapped in its initial angle θ inf. = π/2 − δ by the Hubble friction (see section 7.1). The QCD axion dark matter literature usually employs an axion potential of the form 1 − cos(θ) unlike our choice of cos(θ) from (2.35). The constant piece is irrelevant for dark matter and the correct sign can be obtained by defining θ i = θ inf. ± π, where we take the plus sign for concreteness. However since the Diraxion is light during inflation, its initial angle will receive corrections from quantum fluctuations giving rise to isocurvature perturbations [141] Using H I ≃ 6 × 10 13 GeV and 10 16 GeV < S i < 10 19 GeV one can deduce that 10 −5 < H I /(2π) < 10 −2 , indicating that we can not take the initial misalignment angle to be arbitrarily small. However parametric resonance occurs (see 5.5) for a large part of our parameter space and it has the effect of randomizing the initial misalignment angle, so that the estimate from post-inflationary symmetry breaking with an average angle of π/ √ 3 applies [40]. Since this angle is much larger than the contribution from the fluctuations we take in the following. The Diraxion oscillates around the minimum of its cosine potential (2.32), which reduces to the Diraxion mass m 2 a in the small angle limit. During radiation domination oscillations start at the temperature  which normalized to the entropy density reads [31] ρ mis. a (T a osc. ) s(T a osc. ) Pl.
(5.8) and the relic abundance is found to be Ω mis.
in terms of the the critical density ρ c = 3/(8πG N )H 2 0 , the Hubble rate today H 0 ≡ h · 100 km/s Mpc −1 with h ≃ 0.7 and the entropy density today s 0 related via ρ c /(s 0 h 2 ) ≃ 3.6 eV. The corresponding parameter space was displayed in figure 5. Since we are typically interested in f a = O(10 5 − 10 6 GeV), one can deduce that standard misalignment can not be responsible for the observed dark matter relic abundance.

Topological Defect decay
Cosmic strings are formed via the Kibble mechanism [180][181][182], when the underlying U(1) D symmetry is spontaneously broken down to Z N . Detailed numerical simulations of the cosmic strings decaying to axions e.g. [183][184][185][186][187][188][189][190][191][192][193][194][195][196][197][198][199][200][201] typically have large uncertainties and often disagree with each other, see e.g. [202] for a recent overview. This is why we resort to simple order of magnitude estimates. Cosmic strings are characterized by their string tension or energy by unit length [203] µ str. = π(N f a ) 2 Log m S √ ξH , (5.11) where ξ is a dimensionless length parameter that needs to be determined from numerical simulations. Following reference [190] we take ξ ≃ 1 and note that this parameter could very well be larger. Cosmic strings have a core diameter d ∼ 1/m S and typical distances of L ∼ 1/H, which acts as a regulator for global strings. We obtain the energy density via dividing the string tension by an area, which is given by the typical Hubble volume 4/3π · 1/H 3 over the string separation 1/H to be Apart from cosmic strings there is a second type of defect that might be formed: When the Diraxion oscillates at T a osc. and settles into one of its N degenerate minima the residual Z N is explicitly broken. Consequently domain walls, whose energy density interpolates between the different vacua, will be formed around the time of T a osc. . A similar argument can be made for the case of kinetic misalignment: Since the Diraxion start its oscillation near the top of the potential, even a small fluctuation could lead it to fall into one of the available minima [40]. However these effects typically only occur for scenarios where the global U(1) D symmetry is broken after inflation: The observable universe consists of many patches with randomly distributed, different values of θ i [204] and domain walls separating different patches. Our scenario on the other hand involves (spontaneous and explicit [205]) symmetry breaking during inflation, so one would expect a uniform θ i throughout the observable universe [204]. However this does not guarantee the absence of domain walls. Parametric resonance (see section 5.5) grows fluctuations and the resulting large fluctuations facilitate a non-thermal restoration of the U(1) D symmetry [206][207][208][209][210][211][212] so that the angular field becomes randomized instead of being stuck in its value θ inf. (see the beginning of 5.2). Once the fluctuations are smaller than f a , the U(1) D symmetry is broken again with different initial angles for each Hubble patch leading to domain wall formation at T a osc. [40]. If the condensate was thermalized before parametric resonance can occur, this outcome could be avoided. However as we argue in section 8.1, such an early thermalization is not possible from the Dirac-Weinberg operator, as this UV-dominated rate would dampen the oscillation too much before a rotation could be induced. The only way out is to consider the Seesaw messenger fields as bath particles, which is the scenario sketched in E. Since the symmetry is restored only to be broken again, cosmic strings will also form and the resulting hybrid network of decays can decay for N = 1 [213,214]. A second reason to expect domain walls is that the isocurvature fluctuations can experience a power law growth [23]: Due to the isocurvature fluctuations the different patches that make up our visible universe start with different initial velocitiesθ i and evolve to different field values. By the time T a osc. , when the Diraxion mass becomes cosmologically relevant, domain walls will form to interpolate between the different θ. The surface tension or energy via unit area of a domain wall reads [191] σ DW = 8m a f 2 a (5.13) and one finds the energy density from the ratio of the surface density over the typical length scale. Here we assume O(1) domain walls in our Hubble volume. For a Hubble volume 4/3π · 1/H 3 and a typical surface area of 4π 2 · 1/H 2 the result is (5.14) Domain walls would start to dominate a previously radiation dominated universe at a temperature of The network collapses under its tension when [203] H dec. = µ str. As long as Log (m S /H dec. ) ≳ 8/(3πN 2 ) we find that H dec. ≲ 3m a which means that the axions produced from the hybrid defect decay are produced at a time comparable to the conventional misalignment scenario [203]. For a QCD axion it was found that the logarithm can be as large as 70 [202], so the previous conclusion is likely to hold true. Since the domain walls are expected to form at H ≃ 3m a and decay around the same epoch, they do not have enough time to dominate the energy budget of the universe (compare the formation and domination temperatures in (5.6) and (5.15)). We find that the total energy density of the produced Diraxions is comparable to the conventional misalignment contribution 5.7. The typical energy of Diraxions radiated by the network is about three times their mass [203] , so they are cold dark matter.
For N ̸ = 1 we have more than one domain wall. For odd dimensional operators responsible for the Diraxion mass we sketch in appendix C how to construct such operators while maintaining a domain wall number of one. Since we are however interested in an even dimensional operator with N = 6, we have to invoke an additional operator with a dimension that is co-prime with N = 6, explicitly breaking the residual Z 6 that protects the domain wall [175]. This additional bias term will then lead to the annihilation of domain walls in neighboring vacua [107]. For group theoretical solutions to the domain wall problem see [215,216]. We use the dimension seven operator from (2.35) to define ∆m a which is much smaller than m a from the dimension six operator ∆m a m a ≃ 8.4 × 10 −7 · c 7 c 6 · N f a 10 6 GeV , (5.18) as long as we assume c 6 ≃ c 7 . It is evident that the shift in ε defined in (3.31) due to such values of ∆m a /m a is negligible. The explicit symmetry breaking potential or bias term V 7 induces a different energy for each of the N vacua and the pressure difference between neighboring vacua is proportional to where δ 7 is the phase of the coefficient c 7 and we will take the cosine term in brackets to be of order one from now on. The domain walls collapse if the pressure difference is larger than the domain wall tension implying Since H bias ≪ m a the walls decay long after the onset of coherent oscillations. Reference [175] obtained that the domain walls annihilate before they dominate the energy budget as long as ∆m a m a ≫ 3.3 × 10 −12 · N f a 10 6 GeV , (5.21) which is compatible with (5.18). For our benchmark m a = 0.1 eV and ∆m a /m a from (5.18), we find that the domain walls disappear before the onset of BBN at T ≃ 1 MeV.
There is an upper bound on the bias term from demanding that the vacua on both sides of the domain wall percolate sufficiently, which reads ∆m a /m a < 0.2 [217][218][219] and is always satisfied for our scenario. The energy density of Diraxions is given by (5.14) evaluated at H bias and multiplied by N to take the multiple walls into account. In terms of the energy density one finds [175] ρ bias = 3N ∆m a m a that can be much larger than the previous abundance (5.17) and the misalignment contribution (5.7) due to ∆m a ≪ m a . The numerical simulations of reference [203] revealed that the energy of the Diraxions from such long lived domain walls is about twice their rest mass, so they would be cold dark matter. Fixing the relic abundance would require ∆m a m a ≃ 1.4 × 10 −9 · m a 0.1 eV · f a 10 6 GeV which is typically too small to be realized in our case, see (5.18) with c 6 ≃ c 7 . In other words we expect the Diraxions produced from long-lived domain walls to contribute at most Ω DW a h 2 ≃ O(10 −4 ) due to (5.18). The amount of gravitational radiation produced in the domain wall decay computed by [220] turns out to be negligible for our range of parameters.
Note that string-wall networks with N > 1 can have spherically collapsing closed walls leading to the formation of potentially over-abundant primordial black holes (PBH) [221][222][223]. Reference [224] points out that more numerical simulations are needed in order to conclusively treat this effect and further mentions that the presence of angular momentum could lead to a sufficient deviation from the spherical symmetry impeding the PBH formation. Since here the rotating Diraxion background constitutes a source of angular momentum our scenario could be safe from this effect, but in any case a detailed study would be needed.
We conclude that topological defects do not pose a cosmological problem in our scenario. Since we focus on decay constants in the O(10 5 GeV − 10 6 GeV) range, we find that the contributions from coherent oscillations and topological defect decay are negligible compared to the Kinetic Misalignment and Parametric resonance scenarios encountered in the next sections.

Kinetic Misalignment and fragmentation
If the Diraxion velocityθ is larger than 2m a [31] the rolling pesudoscalar will get trapped in its cosine potential later than for conventional misalignment. As a consequence of its large kinetic energy it does not just probe the harmonic part of its potential (small angles) but instead can start near the hilltop. Due to these effects one obtains the dark matter yield [22] from (3.38) and (3.39) where the factor of two, found numerically in [22] and analytically in [225], encodes the enhancement from the anharmonicity. This mechanism reproduces the observed relic abundance for [31] Ω KM  We can accommodate the parameter region 0.1 eV ≲ m a ≲ O(1 eV) and f a ≃ 10 5 GeV for detecting Diraxion dark matter with PTOLEMY via decays to neutrinos depicted in figure 5 via kinetic misalignment. One can fix m a , f a by adjusting the combination of parameters ε, S i (ε, T RH ) for rotations initiated during radiation domination (reheating). This typically leads to ε ≳ 0.1 once S i is determined via thermalization (see section 8) and T RH is accounted for by (3.10) for oscillations before radiation domination. We picked a rather large value of ε ≃ 0.8 in the above estimate, to ensure the absence of parametric resonance (see (5.39) in the next section), which for previous parameter range would overclose the universe (see (5.41)) with too warm dark matter (consult (5.46)). The used value of m S /(N f a ) is about a factor of three to large to comply with the bound from isocurvature fluctuations for a Hubble induced mass in (7.16). Since all of our estimates come with O(1) uncertainties, we treat this parameter point as right on the edge of the allowed parameter space. If (for rotations during radiation domination) we decrease S i and hence the relic abundance by a factor of three, we can make the Saxion mass compatible again. This parameter region only produces dark radiation (see section 6) with ∆N eff. ≲ 0.01 due to the large m S (see (4.17) and also figure 6) required by the constraint on ε = 0.8.
Alternatively fixing the ratio f s /S i via the relations (4.17) and (4.21) together with explaining the observed baryon asymmetry via Dirac Lepto-Axiogenesis (see (4.18) and (4.22)) allows us to estimate the required Diraxion and Saxion masses for DM from kinetic misalignment:

RH (5.31)
In section 9 we check this parameter range against the constraints from the conditions ε < 1 in (3.33),θ < T from (3.47) and S > T from (3.48).
The Diraxion gets trapped in its potential at the temperature [225] T * ≃ 0.23 GeV · 100 g * (T * ) 1 3 · m a 10 meV defined as the time when the Diraxion's kinetic and potential energy coincide. Kinetic misalignment occurs for T * < T a osc. and one finds that this implies [225] f a < 6 × 10 11 GeV · g * (T * ) 100 · Ω a h 2 0.12 . (5.33) Since the Diraxion scans its potential for a long time, parametric resonance from the Diraxion self-interactions becomes possible leading to a fragmentation of the zero mode condensate into higher momentum excitations. Parametric resonance was first discussed in the context of (p)reheating [226,227]. This effect was then applied to the study of relaxions and axions from monodromy [228][229][230] and recently to kinetic misalignment in [225]. The basic idea is that the the mass term for the higher momentum fluctuations a k is time dependent and acts like the external force for a driven oscillator [32] a k + k 2 + m 2 a cos(θt) a k = 0. (5.34) In the limitθ ≫ m s one finds a narrow resonance band around the momentum k ≃θ/2 with a relative width ∆k/k ≃ m 2 a /θ 2 and that the fluctuations are produced with a rate Γ a PR ≃ m 4 a /θ 3 [230]. It turns out that the Diraxion abundance including fluctuations coincides with the zero mode estimate in (5.25) [35], which can be understood by noting that the characteristic energy scale of both processes isθ [28]. Fragmentation can even lead to a slight enhancement of the relic abundance, because the zero mode redshifting as ρ θ ∼ 1/a 6 is converted into fluctuations that redshift slower than 1/a 6 [225]. As a rule of thumb fragmentation is more important for smaller decay constants [225]: One can show that the condensate does not fragment completely as long as f a < 8.79 × 10 10 GeV [225], meaning that most of its energy density is still in the zero mode and the fragmentation ends after the trapping has taken place. Complete fragmentation occurs before the trapping can take place [225] and implies f a < 4.39 × 10 10 GeV. The analysis of [225] breaks down if the backreaction of the fluctuations on the zero mode occurs before the onset of fragmentation, which can be expressed as m a /H(T * ) < O(10 3 −10 4 ) corresponding to f a < 2.25 × 10 10 GeV. This does not mean that the corresponding parameter space is excluded, just that a more complicated analysis e.g. a lattice study is needed. The aforementioned regions in which kinetic misalignment and fragmentation take place are shown as orange lines in figure 5. Since the typical momenta produced during parametric resonance can be much larger than the Diraxion mass, the fluctuations will be less cold than the zero mode oscillations. We need to ensure that the dark matter matter redshifts enough so that it is sufficiently cold by the time of matter-radiation equality. Reference [225] found that the momentum modes m a a(T * ), with a(T * ) ∼ 1/T * the scale factor at the trapping temperature, undergo the most efficient growth from parametric resonance. The warmness bound on the dark matter velocity reads [231,232] v a eq. ≃ m a · a(T * )/a eq. m a ≲ 2 × 10 −4 . implying T * ≳ 5 × 10 3 T eq. ≃ 5 keV and f a > 9.3 × 10 −2 GeV · g * (T * ) 100 · 10 meV m a · Ω a h 2 0.12 , which is hardly constraining . One should also take into account that, if the field is still not trapped at the time of Big Bang Nucleosynthesis (BBN) then its energy density scaling as ρ θ ∼ 1/a 6 , which is different from ordinary dark radiation ∼ 1/a 4 , could modify the expansion history of the universe and alter the produced light element abundances. This condition T * ≳ 20 keV leads to the bound [225] f a > 8.1 × 10 −2 GeV · 10 meV m a · Ω a h 2 0.12 , (5.38) comparable to the structure formation one.

Parametric Resonance from Saxion oscillations
The non-perturbative process of Parametric resonance [226,227] is analogous to stimulated emission and can also be driven by the Saxion oscillations [55,56]. One finds that the oscillating zero mode of the Saxion field is rapidly converted into higher momentum Saxion and Diraxion fluctuations. Production of QCD axion dark radiation from Saxion parametric resonance was studied in [233,234]. The fluctuations in the angular mode can be produced because both equations of motions in appendix A are coupled as long as S > N f a . Since the energy density in the fluctuations is comparable to the one in the zero mode, these large fluctuations will lead to an effective mass squared for the radial mode that is positive, hence the original U(1) D symmetry gets non-thermally restored [206][207][208][209][210][211][212]. Once the amplitude of the fluctuations redshifts below f a the symmetry is broken again. Saxion fluctuations will later be thermalized and only the rotation together with the Diraxion fluctuations will remain. It was found that the transfer from the condensate into the higher momentum modes begins shortly after the start of the oscillations and for a quartic potential during radiation (matter) domination it ends at around S ≃ 10 −2 (10 −4 )S i [210,212,235], when the backreaction from scattering of the fluctuations with the condensate and themselves becomes important. Parametric resonance from Saxion oscillations requires a violation of the adiabaticity condition |ṁ S (S i )/m S (S i ) 2 | < 1 which can be translated into the condition [23] ε < 0.8 .
which is why the abundance of Diraxion fluctuations is formally equivalent to the result for kinetic misalignment (up to a factor of two) under the replacement ε → 1/2 in the definitions of the yields (3.38) and (3.39). This is also the reason why parametric resonance becomes the leading production channel for dark matter for ε < 1/4. We estimate that   In section 9 we check this parameter range against the constraints from the conditions ε < 1 in (3.33),θ < T from (3.47) and S > T from (3.48).
Owing to the fact that the relativistic fluctuations are produced with momenta k ∼θ(S i ) ∼ m S (S i ) warmness constraints become relevant again: Even though the fluctuations are produced with much larger momenta compared to the fragmentation from Diraxion selfinteractions, the non-perturbative production occurs much earlier at T osc. ≫ T * so in principle the fluctuations can have enough time to redshift sufficiently. The constraint on the velocity is [231,232] v a (T = 1 eV) ≃ k a (T = 1 eV) m a ≲ 2 × 10 −4 (5.44) and signals from the cosmic 21-cm lines are expected to test v a (T = 1 eV) ≳ 10 −5 [236]. We estimate the momentum k a (T = 1 eV) at matter-radiation equality T eq. ≃ 1 eV following [23,40]: Since n a ∼ k 3 a ∼ 1/a 3 one deduces that n a /k 3 a is an adiabatic invariant. From k a ∼ m S (S) and n a ∼ m S (S)S 2 we obtain With this expression we trade the momenta for n a (T = 1 eV) from the yield (5.26) that reproduces the dark matter relic abundance and find One can see that the parameter region for DM from parametric resonance and Dirac-Lepto-Axiogenesis requires larger m S and is thus only viable with early thermalization (see section 8.1).
If we only try to explain the dark matter relic abundance detectable via PTOLEMY (see figure 5), we can eliminate m S via (5.46) for cold enough DM and fix m a , f a together with N = 6 to obtain a lower limit of Here for the rotation during radiation domination we find only one free parameter S i . The above parameter choice corresponds to a Saxion mass of and we had to pick a value of f a > 10 5 GeV, outside of the band potentially detectable with PTOLEMY, to ensure that ε stays below unity (see (3.33)) Cold enough Diraxion dark matter from parametric resonance is not detectable with PTOLEMY via decays to neutrinos. Due to the larger f a we expect ∆N eff. < 0.01 for this range of m S (see section 6 together with (4.17) and the plot 6).

Freeze-In scattering from Dirac Weinberg operator
We estimate the energy density of ν R as follows ρ ν R ≃ ⟨E ν R ⟩ · 2 · 3 · n ν R (6.3) Figure 6. Saxion parameter space in terms of its mass m S and decay constant f a for the isocurvature and dark radiation constraints. The gray region illustrates the amount of dark radiation produced from out-of-equilibrium Saxion decays. The green lines denote the amount of ν R dark radiation produced from scattering with the Dirac Weinberg operator. The dotted (dashed) line corresponds to ∆N eff. ≥ 0.01 (0.1). In the region labelled "Thermalization" the freeze-in approximation breaks down as the ν R would thermalize (see the discussion before and after (6.11)).
where the factor of two counts the helicity states and the three comes from the number of generations. The typical energy of a ν R produced out of thermal equilibrium reads Freeze-in from a thermal bath leads to typical momenta of 2.5 T [254], whereas production from an non-thermal condensate involves momenta of the order of the oscillation frequency m S (S) (see section 5.5). Since we have m S (S) ≪ T (see section 4.1) we take ⟨E ν R ⟩ ≃ 2.5 T . As argued in section 4.3 the Dirac-Weinberg operator will only operate in the Freeze-in regime, which is why we can simplify the Boltzmann equation for the right handed neutrino abundance in terms of the thermally averaged cross sections ⟨σ|v|⟩ and the equilibrium number densities d dt n ν R + 3Hn ν R ≃ ⟨σ|v|⟩ LS→H † ν R n eq. L + ⟨σ|v|⟩ HS→Lν R n eq.
H n S ≃ 7 6 Γ S (S)n eq. L , (6.5) where we made use of the definition of the dissipation coefficient (4.1) and used n eq. H /n eq. L = 4/3. It is straightforward to find the number density frozen-in [244,255]

over a Hubble time 1/H and from that
The scaling relations from 4.3 reveal that the energy density UV-dominated during radiation domination since ρ ν R /ρ SM scales like T . Consequently the production is peaked at the beginning of the oscillations. For production during radiation domination we obtain If the oscillations occur before the completion of reheating, we have to take into account that the entropy of the decaying inflaton heats the SM plasma compared to the decoupled 7 ν R . Another way to come to the same conclusion is the fact that ρ SM is only conserved after reheating. Since ρ ν R ∼ 1/a 4 and the non-relativistic inflaton energy density redshifts as ρ inf. ∼ 1/a 3 we are lead to consider ρ ν R /ρ 4/3 inf. as an invariant, which scales as (1/T 13/3 , 1/T 29/3 ) for (T osc. > T RH > T S , T osc. > T S > T RH ). This implies that production is IR dominated and mostly occurs at T RH where we used the scaling relations to restore Γ S (S i )/H(T osc. ) ≪ 1 and focused on the regime T osc. > T RH > T S where Γ S (S)/H(T ) ∼ T 7/3 , because we typically find that the Saxion reaches its minimum after reheating. The amount of dark radiation is where we eliminated m S using (3.10) to make sure the oscillation starts before the completion of reheating. Comparing equations (6.7) and (6.9) reveals that typically less dark radiation is produced for oscillations during reheating due to the entropy dilution in (6.8).
For both scenarios we find an upper limit of where g * (T osc. ) has to be replaced with g * (T RH ) for production before the end of reheating. Fixing f a via equations (4.17) and (4.21) is equivalent to fixing the amount of dark radiation produced from Freeze-In. The above value for the dark radiation abundance is only a factor of two away from the projected sensitivity of CMB-S4 [250,251], and PICO [252]. CMB-HD [253] would have sufficient sensitivity to probe this scenario. The regions corresponding to observable ∆N eff. from ν R produced via scattering can be found in green in figure 6, where we made use of (4.17). Note that naively extrapolating Γ S (S i )/H(T osc. ) → 1 would lead to ∆N eff. ≃ 0.28, however such large rates would thermalize the ν R (see the region in 6 labelled "Thermalization"). In that case the amount of dark radiation can be found from entropy conservation and it only depends on the decoupling temperature T FO ν R of the three generations of ν R [256] where g ν R = 2 counts the spin degrees of freedom of the right handed neutrinos and the factor of 7/8 arises due to their fermionic nature. As will be explained in section 8, we are not interested in the regime Γ S (S i )/H(T osc. ) → 1, since this would lead to additional damping of the Saxion and a delay in the commencement of its oscillations.
There are additional model dependent processes that can source an abundance of ν R : For the Type I Dirac Seesaw this would be S S ↔ ν R ν R via exchange of a virtual N and for the Type II case one can have L L ↔ ν R ν R by exchange of η. The rates scale as Type II (6.12) and since they are UV-dominated we can ensure that their contribution is negligible by demanding that they are slower than the rate from the Dirac Weinberg-operator in (4.1) evaluated at the largest temperature T max defined in (3.22) and the largest field value S i . By using the Seesaw relations in (2.8) and (2.14) we eliminate the dependence on the heavy messenger masses and find for the Type I Dirac Seesaw that and for the Type II case we arrive at (6.14) which is compatible with the active neutrino mass scale for µ η /Y ν ≳ 6 × 10 15 GeV.

Out-of-equilibrium Saxion decays
The Saxion can decay to Diraxions via its derivative coupling in (2.39) with a width [257] Γ(S → aa) = 1 32π (6.15) and to neutrinos with a decay width, that is given by the expression (5.1) for the decays of Diraxion to neutrinos under the replacement m a → m S . The decay to Diraxions is the dominant mode because m 2 S > 2 ν m 2 ν . Non-thermalized Saxions could basically produce an arbitrarily large ∆N eff. , which is why we require Saxion thermalization (see section 8).
The thermalized Saxions decouple at T D and could decay out of equilibrium to Diraxions at a later time T dec. < T D T dec. = 0.08 Its energy density at T D reads ρ S (T D ) = m S T 3 D ζ(3)/π 2 and since ρ S ∼ n S ∼ 1/a 3 it redshifts to ρ S (T dec. Using this together with the definition of T dec. one finds that [40] ∆N eff. ≃ 0.25 · 10 MeV m S · N f a 10 8 GeV · 100 g * (T D ) · 10.75 g * (T dec. ) 1 12 . (6.17) The current bound ∆N eff. < 0.28 [245] excludes N f a > 10 8 GeV · m S /10 MeV. This bound is not very restrictive to our scenario due to the absence of stellar cooling constraints, which exclude f a < 10 8 GeV for a QCD axion. We showcase this bound as the gray area in figure 6. However the above estimate only holds when the Saxion decouples from the thermal bath while relativistic. In section 8.3 we show that for late thermalization the Saxion can stay in thermal equilibrium until it becomes non-relativistic and hence Boltzmann-suppressed. Furthermore one can see from 6 that the channels for Diraxion and ν R dark radiation exist in different regions of parameter space and only overlap for Saxion masses below the MeV-scale.

Isocurvature perturbations
Since both the Saxion and the Diraxion are light during inflation we expect that quantum fluctuations imprint on them and lead to isocurvature modes [258][259][260][261]. For a quartic potential the typical time-scale over which S relaxes to its minimum N f a is given by t Relax ≃ 1/m S (S) [173,174]. In order for the radial mode to remain stuck in its large initial field value until after inflation, we have to demand that the relaxation time is larger than age of the universe t ≃ 1/H I after inflation with a Hubble rate H I . The corresponding condition immediately implies that there will be isocurvature fluctuations in the Saxion direction, which is typical for Affleck-Dine scenarios.

Dark Matter and Baryon isocurvature
First we investigate dark matter and baryon isocurvature modes. The gauge invariant entropy perturbation reads [261] S a = δ(n a /s) n a /s = δY θ Y θ (7.2) and we compute the amplitude of the isocurvature power spectrum following reference [23] assuming that the Diraxion makes up the whole of dark matter: We find the amplitude of the baryon isocurvature perturbations by an analogous calculation In the above we defined the fluctuations in the angular and radial modes in terms of the Hubble rate during inflation to be Note that in the presence of multiple perturbations, one should add the perturbations in (7.2) coherently to obtain the total amplitude. In our case the baryonic and dark matter fluctuations for N = 6 appear with the same signs, so there is no destructive interference between them, which would give rise to compensated isocurvature perturbations [262]. CMB observations [143] at the pivot scale k * = 0.05 Mpc −1 determined the power spectrum of the adiabatic perturbations to be P ζ (k * ) ≃ 2.2 × 10 −9 and they constrain By rescaling the bound on P a (k * ) to the observed baryon abundance and using h 2 Ω B = 0.0224 as well as h 2 Ω a = 0.12 [245] we find [147] Let us first deal with the contribution from the Diraxion fluctuations: Since only ε depends on the angle θ, we find that the usual and kinetic misalignment scenarios are relevant; parametric resonance depends only the initial Saxion number density. Furthermore we deduce from equations (4.14) and (4.19) that in our scenario the baryon asymmetry always depends on ε. Consequently we find for kinetic misalignment and baryogenesis that where we used (3.31) for the last equality. If we choose θ i + δ = ±π/2 the cotangent vanishes and there will be no perturbations in the Diraxion direction [12]. The Diraxion will be stuck in this value due to the Hubble friction. Hence we have to assume that the Diraxion is not aligned with the minimum of its potential which would instead enforce from equation (3.32) ∂V / D /∂θ ∼ sin(θ i + δ) = 0. For conventional misalignment the amplitude of the dark matter isocurvature spectrum is given by the standard expression [55] with f a replaced by S i P mis.
where ⟨θ 2 i ⟩ was defined in (5.4). Since parametric resonance will randomize the angle we take ⟨θ 2 i ⟩ ≃ π 2 /3 [40] and find One can see that the spectra for parametric resonance correspond to the ones for kinetic misalignment under the replacement N → 4, because here the relic abundance is induced from the Saxion potential V σ ∼ S 4 and not the Diraxion potential V / D ∼ S N [23]. For the baryonic perturbations we find using (4.14) and (4.19) (7.14) We focus on the bound for the dark matter perturbations, since the limit is slightly stronger than for baryons. Our result for the scenario in section 3.1.1, where S i originates from quantum fluctuations, reads 7.15) and the case of a Hubble-dependent mass from the Ricci scalar in section 3.1.2 is constrained to be The corresponding regions were drawn in red in figure 6 and we see that for the case of quantum fluctuations inducing S i , most of the interesting parameter space would be excluded. In section 7.3 we demonstrate under which conditions the Dirac Seesaw models do not induce radiative corrections that violate the aforementioned bounds. The large hierarchy between m S and N f a is the main drawback of using a quartic potential. This problem can be avoided in supersymmetric models such as [40][41][42][43], where the approximately quadratic scalar potential arises from supersymmetry breaking via soft masses in a two field model [263], dimensional transmutation from the RGE running of soft masses [264] or loop corrections in gauge mediated supersymmetry breaking [265] in single field models.

Dark Radiation isocurvature
We define the gauge invariant perturbation in the dark radiation fluid via the relation [266] S DR = 3 4 where we made use of the fact that the photons are adiabatic with respect to the perturbations in the SM plasma to trade ρ γ for ρ SM . The amplitude of the dark radiation isocurvature spectrum reads and we already used the fact that ∆N eff. does not depend on θ in our model, see e.g. (6.7) and (6.9). Dark radiation isocurvature leads to a spatial variation of ∆N eff. compared to its spatial average ∆N eff. [266] ∆N eff. ≃ ∆N eff. 1 + 4 3 P(k * ) DR .

(7.19)
A spatially varying ∆N eff. would lead to spatially varying abundances of light elements produced during BBN. Reference [266] sets constraints on this effect by making use of the 4 He/D data extracted from local galaxies [267] and measurements of D/H obtained from high-redshift Lyman-α absorption systems [268], which correspond to a pivot scale of k * ≃ 1 Mpc −1 and they obtain the limit . (7.20) CMB data with k * ≃ 0.1 Mpc −1 can be used to set constraints of isocurvature modes in the left-chiral SM neutrinos, which can be recast as a bound on dark radiation [266] P DR (k * ) < 10 −10 N eff. ∆N eff.
2 . (7.21) The authors of [269] recast CMB data with variable ∆N eff. , updating a similar study [270] based on WMAP data, and their bound reads eff.

(7.22)
For simplicity we neglect the correlation between the dark matter and dark radiation isocurvature perturbations, which both inherit the fluctuations from S. By setting N eff. ≃ 3 one can see that the bound (7.21) is about an order of magnitude stronger than (7.22), so to be conservative we use (7.21) to set limits on the amplitude of the power spectrum These limits are sub-leading to the ones from dark matter isocurvature in (7.15) and (7.16).

One Loop corrections
Since isorcurvature perturbations typically require a tiny Saxion quartic coupling λ σ where we used the bounds in (7.15) as well as (7.16) and set N − x = 3.5 for the strongest limit, we have to check that radiative corrections to this parameter and to m S are under control. Here we write out only the finite pieces of the loop corrections and neglect logarithmic factors. The mixed quartic coupling with the SM like Higgs, the scalar doublet η for the Type II scenario, or the triplet ∆ in the Type III Seesaw induce corrections of [234] δm where the minus signs take into account that at tree level m 2 S < 0, as well as We find that for both kinds of corrections the mixed quartic couplings need to satisfy |λ σH |, |λ ησ |, |λ ∆σ | < 1.2 × 10 −9 quant. fluct.

Type I-II Dirac Seesaws
For the one-loop correction from the N L − ν R loop in the Type I Seesaw we can recycle the result for a Majorana Seesaw [271,272], because only one chirality of N runs in the loop Both contributions come with a minus sign from the closed fermion loop, which for the two-point function was already absorbed in the definition of m 2 S < 0. The bound from the correction to the negative m 2 S can be re-expressed by using the Type I Seesaw relation in (2.8) as and the one from the quartic is given by It is evident that the bound on M N from the correction to m 2 S is the stronger one and would lead to values of M N that are in conflict with our cosmological assumptions in (3.20) and (3.21), so that the heavy N could be produced from the plasma. On the other hand the bound from the correction to the quartic is compatible with our choices of e.g. T RH = (10 14 − 10 15 )GeV. One way to avoid this conclusion is to assume that there is an accidental cancellation between the corrections δm . This would lead to which is still too small for our purposes. However we can ameliorate this problem by assuming the simultaneous presence of two or more different Dirac Seesaws. The trilinear term connecting σ, η and H in the Type II Dirac Seesaw induces where we eliminate κ by using the Type II Seesaw relation in (2.14) and set a bound on µ η . The limit from the correction to m S is given by and is stronger than the bound from the correction to the quartic (before diagonalization), which is known as a Hybrid-Seesaw [273,274], their one loop corrections to m 2 S could cancel each other as long as (7.37) owing to the fact that the bosonic and fermionic contributions have different signs. In this case the bounds from the quartic in (7.32) and (7.36) also get weakened. If we switch on both types of Dirac Seesaws we find a two loop correction of By itself the Type I (II) Seesaw contribution would give a bound similar to (7.31) ((7.35)) with (1)) and an additional loop factor. Owing to the fact that we have to coherently add the two loop contributions from both Seesaw schemes (there are four ways to glue together the diagrams in figure 1 at two loop) there could be another cancellation: We find that the two loop correction becomes zero if the phases of the different contributions to m ν satisfy δ II = −δ I − π, which together with (7.37) would imply |Y L | ≃ 4|Y ν |.

Type III Dirac Seesaw
Analogously we find for the quartic term connecting σ, ∆ and H for both versions of the Type III Seesaw As was discussed below (7.32) and (7.36), the new degrees of freedom in the Type III scenario even without the aforementioned range of v ∆ are typically so light that we can not avoid their presence in the plasma anyway, which will be exploited in appendix E.

Thermalization
The Saxion only couples to the bath via the Dirac Weinberg-operator, leading to a UVdominated rate (4.1), and the mixed quartic coupling with the Higgs, that leads to a IR-dominated rate to be discussed in the next subsection. We need to ensure that the Saxion is thermalized to avoid the overproduction of dark radiation or relic Saxions.

Early Thermalization
If the Dirac Weinberg-operator is supposed to be fast at early times T ≲ T osc. to thermalize the Saxions, it would be already fast at T osc. . As we saw in section 4.3, an efficient charge transfer from the condensate to the bath would require a rate Γ S (T osc. ) that is a thousand times faster then the Hubble rate for ε = 0.1. Such a large dissipation rate would naively lead to immediate evaporation of the condensate. Reference [157] however showed that since the dissipation rate depends on the oscillating field value (and not the amplitude) and further since it is the coefficient ofṠ in the equation of motion (see appendix A for the coupled equations of motion for both fields) it vanishes twice every period: Once when Γ S ∼ S 2 = 0 and once whenṠ = 0. The authors of [157] found from numerical simulations that such a term Γ S ∼ S 2 does not quickly evaporate the condensate, but instead leads to stronger damping of the amplitude than Hubble friction alone. The backreaction of particle production on an oscillating condensate can be captured by multiplying the amplitude with a factor of exp(−Γ S /H) [275]. A stronger decrease of the amplitude at early times could however decrease the amount of angular rotation produced during the first couple oscillations in (3.30). We also expect that the additional friction from the coupling to the thermal bath would delay the commencement of the oscillations until m S (S) ≃ Γ S (S) H [276], which is below the usual oscillation temperature defined in (3.8). Since the production rate Γ S depends on T , the baryon asymmetry production would consequently be even more inefficient. Studies of oscillating QCD axions with additional friction from a thermal bath [276,277] also find that the oscillations are damped and delayed. We conclude that early thermalization via an effective operator is not necessarily viable and requires a dedicated numerical simulation. Reference [42] considers thermalization of Saxions with Higgsinos via an effective operator similar to (4.1) and avoids the aforementioned problem in the following way: The Higgsino The gray dotted line for the Kaon branching ratio was obtained from [278] and the LHCb constraint continues (with interruptions) until below 5 GeV. Equation (8.19) leads to the purple BBN bound from light Saxion decays to electrons below the muon threshold.
mass gets a correction from the same coupling leading to the thermalization operator and the Higgsinos only become kinematically accessible at some time after the start of the oscillations. We sketch a scenario based on using the potentially light triplet scalar in the Type III Dirac Seesaw (see 2.2.3) in appendix E. Early thermalization before the Saxion has reached N f a has the appealing advantage that the Diraxion gets automatically thermalized as well [40] as a consequence of their coupled equations of motion in appendix A, which removes the warmness-bound for parametric resonance dark matter in (5.46).

Late Thermalization
Here we adopt the thermalization scenario of references [23,40,56], which utilizes the mixed quartic coupling between the singlet scalar and the SM Higgs and summarize the required parameter space. The Saxion can thermalize via the processes SH ↔ HZ, HW and the corresponding rate is [23,40] where we take the weak fine structure constant to be α 2 (T ) ≃ 1/30. Laboratory searches for light scalars constrain the mixing angle θ SH between the Saxion and the Higgs. This angle is defined in (B.4) of the appendix and in the small angle approximation together with m S ≪ m h we can re-express it as Higgs to invisible decays constrain λ σH < 0.01 [279]. We also use the limits from Kaon decays at NA62 [280][281][282], recasts of PS191 data [283] as well as a recast [284] of old CHARM data [285], which are relevant for Saxion masses below about 300 MeV. Additionally we employ the bounds from displaced vertex searches of B-meson decays B → KS(µµ) at LHCb [286,287], where S decays to µµ, which are sensitive below around 5 GeV. Note that these limits assume a 100 % branching fraction of the Saxion to muons and including the kinematically allowed hadronic decay modes above the two pion threshold leads to less stringent bounds [286,287], which is why we treat the aforementioned limits as conservative constraints. For heavier singlets there exists a LEP search [288] for the process e + e − → Z * S, which only excludes θ SH ≳ 0.1. For an overview of the existing laboratory constraints consult [278,289]. The bounds from SN1987A on light scalars were recently reevaluated in [290].
One needs to check that the required values of λ σH do not upset the small tree-level coupling λ σ , whose value is dictated by the isocurvature constraints from section 7.3. Due to one loop corrections to the quartic we require that λ σH < 1.2 × 10 −9 quant. fluct.
1.8 × 10 −4 H-mass from R . (8.4) In the limit m h ≫ m S one would also find that integrating out the SM Higgs would lead to the following tree-level threshold correction to the Saxion quartic [291][292][293] The upper limit in θ SH < 3 × 10 −3 from LHCb [286,287] (see figure 7) can be recast by using (8.3) into (8.6) One can deduce that the corresponding λ σH complies only with the limit for Hubble induced initia field value in (8.4). Furthermore one can see from (7.26) that the tree-level correction is subdominant to the limit on λ σ only for the case of a Hubble induced mass (λ σ < 10 −10 ), but not for the scenario with S i from quantum fluctuations (λ σ < 10 −20 ).

Thermalization during radiation domination
Below T S the Saxion redshifts like non-relativistic matter S = N f a (T /T S ) 3/2 and would lead to an era of matter domination starting at a temperature if it was not thermalized beforehand. During radiation domination we find that Γ S,H /H ∼ 1/T so thermalization is IR dominated. The thermalization temperature is found to be (8.9) To avoid complications from the oscillatory effective masses of the Higgs for S ≫ N f a we work in the limit T th < T S [23] with T S defined in (3.12), which implies The condition T th > T M on the other hand leads to . (8.11) One should keep in mind that the SM like Higgs has to be as abundant as radiation for this process to work, which is why we require T th > m h = 125 GeV θ SH > 1.4 × 10 −7 · g * (T th ) 100 1 4 .

(8.12)
The Higgs could also receive a mass correction from its coupling to S namely ∆m 2 H ≃ 2λ σH N f a S, where we dropped the S 2 term following [23,40], because it is expected to be subdominant for T < T S owing to S = N f a (T /T S ) 3 2 < N f a . Demanding that T th > ∆m H gives an upper bound Additionally we need to ensure that the Saxion does not receive a large thermal correction from its coupling to the abundant Higgses m S (S i ) > √ λ σH T osc. which can be cast as 14) The bounds in (8.10) and (8.11) are the strongest thermalization constraints and we depict them in figure 7 by fixing f a /S i via (4.17), which explains the dependence on m 2 ν and ∆N eff. . We find that all the available parameter space would be excluded by the constraints from meson decays unless we take S i ≲ 0.1M Pl. . Furthermore for this initial field value we have to require that ∆N eff. > 2.8 × 10 −3 , or else thermalization before Saxion domination would be excluded by LHCb. The original Lepto-Axiogenesis parameter space for a quartic potential involves Saxion masses that typically lie below the GeV-scale, whereas we will show in (9) that our Saxion can be heavier.

Thermalization during reheating
For the case where the oscillations begin before the completion of reheating our analysis finds that typically T S < T RH which together with T th < T S leads to the conclusion that thermalization will again proceed during radiation domination. Stil we need to demand For our parameter space this constraint is subdominant compared to (8.10) and (8.11).

Thermalized Saxion decays
In section 6.2 we showed that a significant part of the parameter space would be excluded by dark radiation constraints if the thermalized Saxions decayed long after their decoupling from the bath. Hence we require that the Saxions are in thermal equilibrium when they become non-relativistic. The mixing with the SM like Higgs allows for the following decay modes to the light SM leptons with l = e, µ, (8.16) where we neglect the phase space suppression for our first estimate and the corresponding decay temperature reads BBN will not be affected if the Saxions decay before the neutrino decoupling at T = 2 MeV implying This limit was depicted in figure 7. Additionally the Saxions need to be Boltzmannsuppressed when the SM neutrinos decouple, because else at T < T dec the Saxions are still kept in equilibrium via decays and inverse decays to electrons and inject entropy into the SM plasma. This process would heat only the SM bath so the already decoupled neutrinos are cooled leading to a reduction in ∆N eff [294]. Reference [23] obtained that the lower bound ∆N eff < −0.44 [295] is only compatible with (see also [296] for a similar analysis involving heavy QCD axions) As it turns out our scenario typically requires Saxions with masses around the GeV-scale (see sections 5.4 and 5.5), which means that on the one hand decays to mesons are possible and on the other that the Saxion would be Boltzmann-suppressed at BBN. The decay to SM fermions ψ is the dominant mode compared to Diraxion final states as long as From the isocurvature limits (7.15) and (7.16) one can deduce that the right hand side in the above is of O(10 −5 − 10 −11 ), which is much smaller than the coupling to SM states for GeV-scale Saxions.

Maximum Yield
The thermalization temperature allows us to determine an upper limit on the charge yield [22,40,42,43] where we used T th. < T S so that the result only depends on m S . This relation can be understood as follows: During thermalization the energy density of the oscillations ρ S = (1 − ε)m S n S is dumped into the bath [37] and only the pure rotation with ρ θ = εm S n S will remain. The entropy density released during thermalization is found from the first law of thermodynamics to be s f ∼ ρ S /T th and consequently we obtain for the dimensionless dilution factor 8 that leads to Y θ /∆ ∼ T th /m S . The entropy generation is maximal for the case of Saxion domination, for which an equal sign in (8.22) would apply. We express the maximum charge yield in terms of the mixing angle θ SH < 3 × 10 −3 by using (8.9) together with (8.22) and ε ≃ 1 (see the discussion in the beginning of 5.4 and the next section). On the upper axis of figure 6 one can read off the value of Y max θ for a given m S .

Discussion
Late Saxion thermalization via the Higgs portal before Saxion domination requires S i ≲ 0.1M Pl. (see the plot 7) and for concreteness we saturate this value. This leads to a * ★ Figure 8. Parameter space for the cogenesis of the baryon asymmetry and dark matter relic abundance via kinetic misalignment or parametric resonance spanned by the Diraxion mass m a and the Saxion mass m S . The decay constant f a was determined via (4.17) in terms of ∆N eff. and m 2 ν . In the upper plot we fixed ∆N eff. = 1.4×10 −3 and the point marked with a star corresponds to (m S , m a , f a ) = (500 GeV, 1 eV, 1.1 × 10 6 GeV). Here the Diraxion has a lifetime of 4 × 10 4 in units of the age of our universe. In the lower plot we fixed ∆N eff. = 5 × 10 −3 and the point marked with a star corresponds to (m S , m a , f a ) = (2.2 GeV, 31 meV, 2.8 × 10 6 GeV). * ★ Figure 9. Parameter space for the cogenesis of the baryon asymmetry and dark matter relic abundance via kinetic misalignment or parametric resonance spanned by the Diraxion mass m a and the Saxion mass m S . The decay constant f a was determined via (4.21) in terms of ∆N eff.
(the upper limit during radiation domination defined in (6.10)) and m 2 ν . In the upper plot we fixed the upper limit ∆N eff. high oscillation temperature for the Saxion in (3.8) and (3.14), which is why we typically need reheating temperatures above around 10 14−15 GeV, which is close to the limit T RH < 6.5 × 10 15 GeV inferred from H I < 6 × 10 13 GeV [143] and could be seen as a disadvantage of our scenario compared to Leptogenesis or Lepto-Axiogenesis [40], that both work for reheating temperatures as low as 10 9 GeV [297]. Thermalization with S i = 0.1M Pl. is only viable as long as ∆N eff. > 2.8×10 −3 for Saxions around the GeV-scale. To ensure that the processes encoded in Γ S never thermalize, we fix f a /S i via equations (4.17) (radiation domination) and (4.21) (during reheating), which is equivalent to fixing the amount of ν R dark radiation produced in (6.10). We find that we typically need Γ S (S i )/H(T osc ) < 0.1 to avoid ε > 1 (see (3.26) and (3.33)), so the upper limit on ∆N eff. is never saturated and our cogenesis setup does not lead to observable amounts of dark radiation. For N = 5 the bounds from the self-consistency criteria in (3.47) and (3.48) would be the most relaxed, however we find that cogenesis of the baryon asymmetry and dark matter would take place in regions with m a = O(keV). While we can generate the correct relic abundance in this regime, dark matter will not be long-lived enough due to its decay to neutrinos (see (5.3) for the lifetime). Hence we fix N = 6 for the rest of our analysis, which leads to stable enough Diraxions with m a = O(meV − eV). The downside of this parameter range of N and S i is that we find ourselves in a region, where the warmness constraint (5.46) on parametric resonance dark matter will always be violated. Thus we either have to invoke a scenario for early thermalization of the coupled Saxion-Diraxion system as sketched for the Type III Dirac Seesaw in E or tune the eccentricity parameter ε defined in (3.31) close to ε ≃ 0.8 (consult section 5.5 for details).
In figures 8 and 9 we plot the baryon asymmetry as red line and the dark matter relic abundance from kinetic misalignment as a blue line. Both quantities are overproduced above their respective lines. Parametric resonance dark matter is possible in the blue shaded region surrounded by the blue dashed line indicating Ω PR a h 2 ≳ 0.12 and the dotted line for ε = 0.8. Inbetween these lines dark matter from parametric resonance is overproduced. Our parameter space for oscillations during radiation domination (reheating) can be seen in figure 8 (9). In the gray shaded area parametric resonance dark matter would be cold enough and in the orange regions we begin to violate the self-consistency criteriaθ < T from (3.47) and S > T from (3.48). The region in which the upper limit on the yield from thermalization in (8.22) would not be able to explain the dark matter relic abundance via kinetic misalignment is completely within the orange area. The purple region is excluded by neutrino cooling from too light Saxions during BBN, see (8.20). In the cyan area of figure 8 (9) the Saxion oscillation would take place during reheating (radiation domination) as described by equation (3.10). The first panel in 8 shows a point (m S , m a , f a ) = (500 GeV, 1 eV, 1.1 × 10 6 GeV), which would work with early thermalized parametric resonance (here the reheating temperature would actually be to large compared to the bound inferred from H I ). While m a is of the right order of magnitude to be detected by PTOLEMY [176,177] via decays to neutrinos [178,179] its lifetime of about about 4×10 4 times the age of the universe is far too large, so there would be no appreciable number of produced neutrinos. Note that here we can have ∆N eff. < 2.8 × 10 −3 because the Saxion mass is around 100 GeV and thus not constrained by B-meson decays. The second panel in the same figure shows a point (m S , m a , f a ) = (2.2 GeV, 31 meV, 2.8 × 10 6 GeV), where cogenesis with kinetic misalignment is possible for ε ≃ 0.8. If we try to generate a signal from kinetic misalignment Diraxions decaying to neutrinos detectable with PTOLEMY, we have to choose the same value ε (see the discussion in section 5.4). While ε = 0.8 is allowed by the limits in (3.26) and (3.33), here the S i -dependent Diraxion mass m a (S i ) starts to become comparable in magnitude to the initial Saxion mass m S S(S i ), which might potentially hinder the onset of the coherent rotation. Reference [37] solved the equations of motion for N = 6 and ε = 0.8 numerically and found that a rotation can be initiated, however their analysis only covered an approximately quadratic potential.
This conclusion can be somewhat ameliorated for rotations during reheating as depicted in 9: Again the upper panel showcases a point (m S , m a , f a ) = (11.2 GeV, 50 meV, 2.2 × 10 6 GeV) that would work with early thermalization for parametric resonance dark matter. This panel also demonstrates that kinetic misalignment and parametric resonance scale differently during reheating, which can also be deduced by comparing the second line of (5.26) with the second line of (5.41). One finds that the relic abundance for parametric resonance scales with f a /S i and becomes independent of S i after one uses (4.21). The different scaling allows us to find the point (m S , m a , f a ) = (63 MeV, 1.5 meV, 6 × 10 7 GeV) depicted in the lower panel of 9, where the lines for the baryon asymmetry and kinetic misalignment meet for ε ≃ 0.4. The required Saxion mass of 63 MeV is however excluded, if thermalization is supposed to occur via the Higgs portal, see figure 7. Note that the values of ∆N eff. mentioned in the caption of 9 where computed using the estimate for the maximum amount of dark radiation during radiation domination in (6.10).
For our parameter space of interest we find that f a ≃ (10 5 − 10 6 ) GeV, which means that regular misalignment and topological defect decay can at most only contribute a tiny fraction of the DM relic abundance. Inspection of 5 reveals that for these values of f a and kinetic misalignment one expects a completely fragmented Diraxion, so in principle a lattice study is needed. Consequently the Diraxion will not be a zero mode condensate for either production from parametric resonance or kinetic misalignment. As argued in the beginning of this section, we find ∆N eff. < 0.028 in the regions, where we can explain both the observed baryon asymmetry and dark matter. A detection of a lager amount of dark radiation by next generation experiments would exclude only the cogenesis scenario. Thus such a detection would imply that this construction can be responsible for either baryogenesis or the origin of dark matter.

Conclusion
• The Diraxion from the Dirac Seesaw: We showed that all three versions of the Dirac Seesaw mechanism automatically provide us with a PNGB that we call the Diraxion, whose cosmological implications have previously not been analyzed. This Goldstone mode originates from a global U(1) D symmetry that forbids couplings of ν R to the SM Higgs and we assume that this symmetry is explicitly broken by quantum gravity making the Diraxion massive. Furthermore we impose a separate (gauged or residual) symmetry like e.g. U(1) B-L to forbid all Majorana masses. The Diraxion has no direct coupling to SM fields and its connection to charged fermions or photons arises only at one or two loops respectively, which is why our set-up is unconstrained by stellar and supernova cooling arguments.
• Diraxion rotation: The Diraxion is accompanied by a radial mode called the Saxion, whose vev can undergo a large excursion during inflation and which oscillates in the early universe. We discuss two ways of inducing this large vev with a Hubble-dependent mass during inflation being the least excluded (see figure 6). If the Diraxion mass comes from a higher dimensional operator, it can convert a part of the Saxion's oscillatory motion into a coherent angular rotation around the bottom of the scalar potential. For our parameter space we find that we need a dimension six operator to do so. If this operator involves additional insertions of the U(1) B-L -breaking scalar, we can also generate the cosmologically required mass scale m a = O(meV−eV) without assuming a small Wilson coefficient. Conventionally we express the abundance of Diraxions εn S in terms of the Saxion number density n S and a dimensionless eccentricity parameter 0 < ε < 1. We have to keep ε below unity, because else the Diraxion would start oscillating before the Saxion leading to it being trapped in a false vacuum, which hinders the onset of the coherent rotation.
• Dirac-Leptogenesis: We assume that all of the Seesaw messenger fields are too heavy to be present in the plasma. The Diraxion rotation is thermodynamically stable and can act as a background field enabling spontaneous baryogenesis via Dirac-Leptogenesis [48] from the Dirac-Weinberg operator. Scattering processes like HS → Lν R , LS → H † ν R conserve the total B-L, but not the individual SM and right handed neutrino lepton numbers. This effectively leads to equal and opposite chemical potentials between the SM leptons and the right handed neutrinos µ SM B-L = µ ν R . Since the electroweak sphaleron process only couples to the left-chiral lepton doublet one can convert the asymmetry stored in the lepton doublet into a baryon asymmetry. The leptonic asymmetry is not washed out, because we ensure that the only interaction producing ν R is always out of equilibrium.
• Dark Matter: The Diraxion rotation can be responsible for dark matter via the kinetic misalignment mechanism. Additionally fluctuations of the Diraxion can be produced via parametric resonance from the Saxion oscillation. There is no domain wall problem in this construction, because domain walls can immediately decay to Diraxions either via a string-wall network with a single domain wall or for more domain walls via the well known bias-term from the Diraxion mass. Both contributions from topological defect decay are comparable to the relic abundance from standard misalignment, which is negligible for our range of decay constants f a ≃ (10 5 − 10 6 ) GeV. Kinetic misalignment can lead to Diraxion dark matter with 0.1 eV ≲ m a ≲ O(1 eV) and f a ≃ 10 5 GeV (see figure 5), which can be potentially tested via its decay to neutrinos with future cosmic neutrino background experiments like PTOLEMY. For these parameters we expect that the Diraxion zero-mode condensate is completely fragmentated into higher momentum excitations and in principle a lattice study is needed to treat this regime. Additionally the resulting Saxion mass is in tension with isocurvature constraints, which is why for this range of parameters we can only accomodate a fraction of the relic abundance. Dark matter from Saxion induced parametric resonance is constrained by its warmness. This channel can not produce cold enough Diraxions detectable with PTOLEMY, because we require f a ≫ 10 5 GeV to ensure ε < 1. In order to set the relic abundance via kinetic misalignment and to avoid too much and too warm parametric resonance dark matter we have to set ε ≃ 0.8, which shuts off Saxion induced parametric resonance [23]. This implies that the initial (evaluated at S i ) Diraxion mass is comparable in size to the initial Saxion mass, which brings our scenario closer to the mechanisms relying on heavy axion oscillations like e.g. [11,50]. For a quadratic potential with N = 6 and ε = 0.8 it was shown numerically by the authors of [37] that a rotation can be initiated.
• Parameter space and Cogenesis: Our parameter space is spanned by the Diraxion and Saxion masses m a , m S together with the initital Saxion field value S i and the Saxion vev today N f a , where N is the domain wall number. Thermalization fixes S i ≲ 0.1M Pl. and we can eliminate f a by demanding that the Saxion does not have fast interactions with the bath, which also fixes the amount of right handed neutrino dark radiation. Throughout this work we use N = 6 to keep the Diraxion light and long-lived enough. These inputs allow us to determine the values of m a , m S for generating the correct dark matter abundance together with the right baryon asymmetry. Our parameter space for cogenesis with m a = O(meV − eV) and m S = O(100 MeV − 100 GeV) suffers from too warm Diraxions from parametric resonance, which is why we either need thermalization from additional degrees of freedom in the bath, present e.g. in the Type III Dirac Seesaw (see appendix E), or a region where ε ≃ 0.4 (0.8) for oscillations during reheating (radiation domination) depicted in figure 9 (8). Our Saxion is typically predicted to be heavier than for the standard Lepto-Axiognesis scenario. It could be produced by (heavy) meson decays as depicted in figure 7. Due to the size of S i and m S the Saxion oscillations occur at temperatures around 10 15 GeV, which implies comparable reheating temperatures close to the limit T RH < 6.5 × 10 15 GeV inferred from H I < 6 × 10 13 GeV [143]. This is a drawback compared to regular Leptogenesis or Lepto-Axiogenesis which work for reheating temperatures as low as 10 9 GeV [40,297].
• Dark Radiation: This setup can produce ν R dark radiation with 2.8 × 10 −3 ≤ ∆N eff. ≤ 0.028, where the lower limit applies for S i = 0.1 M Pl. and Saxion masses below around 5 GeV due to thermalization via the Higgs portal. While explaining both the baryon asymmetry and dark matter relic abundance involves smaller values of ∆N eff. , fixing only one of these observables allows us to generate more dark radiation and and to saturate the previous upper limit (see equations (5.30)-(5.31) and (5.42)-(5.43)). There is no observable dark radiation for cogenesis, because this would require ε > 1. Consequently we can use next generation measurements of ∆N eff. to test which cosmological history is realized in this framework.
• Isocurvature constraints: Dark matter and baryon isocurvature constraints enforce a value of m S that is very small compared to the U(1) D breaking scale N f a (see figure 6). For S i generated from quantum fluctuations we need m S /(N f a ) < O(10 −10 ) and for a Hubble-dependent mass we find m S /(N f a ) < O(10 −5 ). We showed that one-loop corrections in the various Dirac Seesaw models do not upset this tuning. This comes at the price that the heavy messenger fields responsible for the Dirac-Weinberg-operator could be potentially light enough to be present in the plasma, so the description in terms of the Dirac-Weinberg-operator breaks down. To avoid this we assume an accidental cancellation between the one loop corrections to m 2 S from the heavy Dirac fermion N for the Type I Seesaw and the heavy doublet η for the Type II Seesaw (see (7.37)). Thus our scenario is a Hybrid-Seesaw [273,274] and and each contribution could be responsible for sourcing one of the two mass splittings observed in neutrino oscillation experiments implying that only one generation of N would be needed. The messenger fields for the two versions of the Type III Dirac Seesaw are usually much lighter than N, η so a separate study is needed to investigate their cosmological evolution and impact on the presented mechanism. Our work predicts dark radiation isocurvature correlated with the dark matter and baryon isocurvature fluctuations, potentially detectable as neutrino isocurvature modes in the CMB. The dark radiation isocurvature signal is subleading if we saturate the bound from dark matter isocurvature.

• Extensions and Outlook:
It would be worthwile to consider (a) the case of thermalized, lighter messenger fields for all Dirac Seesaws, (b) a supersymmetric set-up to ameliorate the strong isocurvature bounds on the non-supersymmetric quartic potential and (c) whether the Saxion could play the role of the inflaton, so Saxion thermalization is automatically obtained from successful reheating. Since each of these aspects constitutes a substantial modification of our analysis, we leave them for future investigation.
the local string is determined from the minimization of the system's kinetic energy leading to 9 [301,302] The important point is that the system of domain walls and two types of strings will be unstable if either of the domain wall numbers is one [301,302]. Since we need n σ > 5 for Axiogenesis [40] this leaves the two choices n φ = 1, ω (φ) σ = 0 or n φ ̸ = 1, ω (φ) σ = (n φ − 1)/n σ . A non-trivial winding number ω (φ) σ implies mixed anomalies between U(1) D and U(1) B-L . Since the gauge symmetry is assumed to be abelian, this does not lead to non-perturbative contributions from instantons to the Diraxion mass. On the other hand, if U(1) B-L is embedded into a larger non-abelian group such as the Pati-Salam hypercolor SU(4) c [303], these unwanted effects reappear and might be even compounded by small size UV-instantons [304,305], which could limit the possible UV-completions of our setup. The dynamics of φ in the early universe can be neglected as long as |λ σφ | S i > Max T RH , T max , H I 2π . (C.5) As long as λ σφ < 0 there is no danger of restoring the B-L symmetry via the large vev S i .

D Initial field value from a coupling to the Inflaton
A second way to generate a Hubble dependent mass contains one (or more) coupling(s) to the inflaton field χ that could take the following forms V (σ) ⊃ c 1 |χ| 2m |σ| 2n that were suggested by references [173,174,306], [307] and [308] respectively. One can understand the effect of these couplings by focusing e.g. on the second term and using the Friedmann equation to eliminate V (χ) for slow roll inflation, which implies an effective mass squared 3c 2 H 2 I . In general the idea is, that the large initial field value χ i of the slowly rolling inflaton field sources an effective mass term that is initially larger than the Hubble rate m S (χ i ) > H I . As the inflaton field decreases the effective mass will become smaller than the Hubble rate after N last e-folds of inflation (counted from the end of inflation) so that quantum fluctuations can push S to the value (we assume m S < H I ) [173,174,306] That implies that isocurvature fluctuations will not be produced before N last . Observations by the Planck and WMAP collaborations only constrain isocurvature modes whose 9 Here we assumed that vσ wraps once around the global string and vφ winds once around the local string momenta are below the pivot scale k * = 0.1 Mpc −1 [143]. The matter power spectrum extracted from Lyman-α data is only sensitive to perturbations at scales of 0.2 MPc −1 ≲ k ≲ 10 Mpc −1 [309]. If the Saxion fluctuation has comoving momenta k ≥ 10 Mpc −1 , current experiments do not set a limit on the isocurvature power spectrum. Reference [306]  where we suppressed subleading terms depending on the temperature today and number of relativistic degrees of freedom. Determining the evolution of the effective mass from the evolution of χ and checking whether N last is realizable, requires specifying an inflationary model and a dedicated analysis beyond this work. It is important to note that we chose only effective couplings between the inflaton and the singlet scalar in (D.1) in order to not destabilize the flat inflaton potential too much [307] and to avoid additional, potentially large, contributions to the energy density of the universe [307]. The second and third operator in the above are motivated by scenarios where the inflaton is protected by a shift symmetry [310], that is only explicitly broken by V (χ). Consequently these operators do not introduce a new source of shift symmetry breaking.

E Early Thermalization from Type III Dirac Seesaw
To realize early thermalization, that occurs some time after the initiation of oscillations, but still early enough to thermalize the parametric resonance Diraxions, we need a particle with a thermal abundance and a renormalizable coupling to S, such as the iso-triplet ∆ from the Type III scenario 2.2.3, which may be much lighter than the other new degrees of freedom (see (2.24)). Early thermalization could proceed as follows: ∆ is relativistic and rapidly scatters with the bath due to its gauge interactions. For this to stay true we need to ensure that the corrections to the triplet mass squared λ H∆ S 2 i (tree level) and λ 2 4 /(16π 2 )S 2 i (one-loop) stay below T 2 osc. , which has typical values of O(10 15 GeV) for successful cogenesis. This implies that λ H∆ < 10 −6 · 0.1M Pl S i 2 · T osc. 10 15 GeV 2 , (E.1) In general the Saxion may also receive a thermal mass from its coupling to ∆. The thermal mass is smaller than the tree level mass as long as Max(λ H∆ , λ 4 ) < S i 0.1M Pl · 10 15 GeV T osc. We expect interaction rate between triplets and Saxions to scale as 10 e.g. Γ ∼ Max[λ 2 H∆ , λ 2 4 ] S 2 /T [153,154] so that during radiation domination Γ/H ∼ S 2 /T 3 ∼ 1/T , which is IR dominated. Hence the scattering with ∆ can be taken to be slow when the oscillations start implying T th. < T osc. and further Max(λ H∆ , λ 4 ) < 1.4 × 10 −3 · 0.1M Pl S i · T osc. 10 15 GeV 2 . (E.4) The thermalization rate peaks while the ∆ are relativistic and then rapidly decreases once they are Boltzmann-suppressed. In order to avoid warmness bounds on the Diraxion abundance from parametric resonance, thermalization during radiation domination needs to occur before the Saxion field value has reached S ≃ 10 −2 S i [210,212,235], hence we require T th. > 0.01 T osc. , from which Max(λ H∆ , λ 4 ) > 1.4 × 10 −4 · 0.1M Pl S i · T osc. 10 15 GeV 2 (E. 5) follows. Comparing this to (E.1) and (E.2) reveals that only λ 4 is large enough for thermalization, which is not in conflict the isocurvature bounds in (7.40) and (7.41) (at least for the Hubble induced Saxion field value). Since ∆ is charged under the U(1) D symmetry, its interactions with the condensate can transmit a fraction of the Noether charge into an asymmetry of ∆. Reference [33] considers the charge transfer between rotating condensates in detail. The asymmetry in ∆ is then converted via its Yukawa couplings (essentially the same reactions as in 4.1 but with S replaced by ∆) into an asymmetry of the left chiral leptons that gets reprocessed into a baryon asymmetry via the sphalerons. Once ∆ becomes Boltzmann suppressed we can integrate it out and recover the Weinberg-operator.