Maximal temperature of strongly-coupled dark sectors

Taking axion inflation as an example, we estimate the maximal temperature ($T_{\rm max}^{ }$) that can be reached in the post-inflationary universe, as a function of the confinement scale of a non-Abelian dark sector ($\Lambda_{\rm IR}^{ }$). Below a certain threshold $\Lambda_{\rm IR}^{ }<\Lambda_{\rm 0}^{ } \sim 2\times 10^{-8}_{ } m_{\rm pl}^{ }$, the system heats up to $T_{\rm max}^{ } \sim \Lambda_{\rm 0}^{ }>T_{\rm c}^{ }$, and a first-order thermal phase transition takes place. On the other hand, if $\Lambda_{\rm IR}^{ }>\Lambda_{\rm 0}^{ }$, then $T_{\rm max}^{ } \sim \Lambda_{\rm IR}^{ }<T_{\rm c}^{ }$: very high temperatures can be reached, but there is no phase transition. If the inflaton thermalizes during heating-up (which we find to be unlikely), or if the plasma includes light degrees of freedom, then heat capacity and entropy density are larger, and $T_{\rm max}^{ }$ is lowered towards $\Lambda_{\rm 0}^{ }$. The heating-up dynamics generates a gravitational wave background. Its contribution to $N^{ }_{\rm eff}$ at GHz frequencies, the presence of a monotonic $\sim f_{\rm 0}^3$ shape at $(10^{-4}_{ } - 10^2_{ })\,$Hz frequencies, and the frequency domain of peaked features that may originate via first-order phase transitions, are discussed.


Introduction
As the nature of dark matter remains unresolved and non-standard ideas have become an accepted part of the speculation, one of the avenues is to envisage the existence of a whole dark sector. There is a great variety of possibilities for the field content of the dark sector and for its interactions with the visible one. Yet any dark sector surely couples to gravity, and then it is natural to assume that it connects to inflationary dynamics as well.
If the dark sector consists of a non-Abelian Yang-Mills theory, so that gauge invariant operators have dimension 4, then its interactions with the Standard Model can be very weak, possibly even suppressed by the Planck mass squared (cf., e.g., refs. [1][2][3] and references therein). Weak interactions between the dark and visible sectors allow the dark sector temperature to differ from the Standard Model one. We would like to know how high the dark sector temperature can be, as this affects several phenomena, such as the spectrum of gravitational waves that gets generated; the efficiency with which dark matter candidates can be produced; and the kind of thermal phase transitions that can be encountered.
When we talk about thermodynamic notions, it is a relevant question under which conditions the temperature can be defined at all. For a non-Abelian Yang-Mills theory, the generic equilibration rate, originating from kinematically unconstrained 2 → 2 scatterings, is of order Γ g ∼ α 2 T dark , where α ≡ g 2 /(4π) is the gauge coupling. An upper bound on the temperature is obtained by comparing this with the Hubble rate of a radiation-dominated expanding universe, H ∼ max{T 2 dark , T 2 visible }/m pl , where m pl is the Planck mass. For T dark < α 2 m pl , the equilibration rate exceeds the Hubble rate, i.e. Γ g > H. If we consider dark sectors with α ∼ 0.3, it is therefore in principle meaningful to discuss temperatures almost up to the Planck scale. On the other hand, if the Hubble rate is dominated by a temperatureindependent part, like a vacuum energy density, there is also a lower bound on T dark . We will return to a posteriori comparisons of the equilibration and Hubble rates.
In order to carry out a concrete discussion, we adopt a specific inflationary scenario that can indeed be argued to thermalize efficiently (cf. sec. 2.7), namely that of non-Abelian axion-like (or natural) inflation [4]. The parameters of the inflaton potential are fixed from standard "cold inflation" predictions, to match Planck data [5]. The heating-up dynamics is characterized by the gauge coupling α that does not affect inflationary predictions at leading order. Apart from the self-interactions of the Yang-Mills plasma, α also parametrizes the interactions between the dark sector and the inflaton, and for this we adopt the form of the pseudoscalar operator, where ϕ is the inflaton field, F c µν is the Yang-Mills field strength, c is a colour index, and f a is the axion decay constant. The advantage of this interaction term is that concrete (even if so far incomplete) information is available about the friction and mass corrections that it leads to. We normally reparametrize α through a dark confinement scale Λ IR , cf. eq. (2.21). Furthermore, for simplicity, we denote T ≡ T dark in the following, and assume that the effect of the visible sector can be neglected in the period of time that we are interested in. 1 Our presentation is organized as follows. After an exposition of our general setup (cf. sec. 2), we first introduce the concept of a stationary temperature. The latter permits for a simple qualitative estimate of the energy density that the Yang-Mills plasma obtains during inflation (cf. sec. 3.1). For a quantitative understanding, we then proceed to a numerical solution of the maximal temperature, which is somewhat higher than the stationary one (cf. sec. 3.2). After elaborating upon physical implications for gravitational waves (cf. sec. 4), we turn to a summary and outlook (cf. sec. 5). 1 If the dark sector consists of a relativistic plasma, this assumption is roughly equivalent to T dark > T visible .
There are concrete scenaria where it has been argued that the dark sector indeed heats up first and injects subsequently a part of its energy density into the visible one, so that the dark sector could be hotter than the Standard Model, see e.g. refs. [6][7][8] and references therein.

Evolution equations
Given the fast thermalization rate of non-Abelian gauge theory (the arguments for this are revisited in sec. 2.7), we carry out our discussion assuming that the notion of a local temperature-like quantity can be defined. The degrees of freedom are then the average inflaton field,φ, and the plasma temperature, T . The equations governing their evolution can be written as (a justification from energy conservation follows below eq. (2.4)) where e r and p r denote the energy density and pressure of radiation. Furthermore H ≡ȧ/a is the Hubble rate; V is the inflaton potential; 2 V x ≡ ∂ x V ; and Υ is a friction coefficient, which transfers energy from the inflaton to radiation degrees of freedom (cf. sec. 2.4). If we set T → 0 and Υ → 0 as initial conditions, the plasma remains at zero temperature, and we return back to normal cold inflation. Denoting the pressure and energy density carried by the inflaton by and multiplying eq. (2.1) byφ, the evolution equation forφ can equivalently be expressed aṡ Summing together eqs. (2.2) and (2.4) yields the overall energy conservation equation, e + 3H(e + p) = 0, where e ≡ e ϕ + e r and p ≡ p ϕ + p r . In the same notation the Friedmann equation reads H 2 = 8πe/(3m 2 pl ), where m pl = 1.22091 × 10 19 GeV.
As we will see in sec. 3.2, the system just defined can cross a first-order phase transition. In this case, eq. (2.2) needs to be supplemented by another equation, valid when the system is in a mixed phase. The technical reason is that in a mixed phase, the temperature stays constant at T = T c , so thatṪ = 0. At the same time, the energy density has a discontinuity, e r (T + c ) − e r (T − c ) > 0, so that ∂ T e r | T =Tc diverges. Therefore a naive evaluation of the time derivative is ambiguous,ė r =Ṫ ∂ T e r = "0 × ∞".
In the real world, a mixed phase can incorporate complicated physics (bubble nucleations, sound wave dynamics, turbulence). However, the overall picture should be well captured by an adiabatic approximation. In this treatment, we re-parametrize e r (t)| T =Tc through a volume fraction, u, as In contrast, the pressure p r is continuous at T = T c , since it equals minus the free energy density, and therefore independent of u. Thereby eq. (2.2) gets replaced witḣ We note that the potential V and its derivatives are well-defined, since the inflaton field does not undergo any phase transition in our setup. The solution of the differential equations needs now to be complemented by a monitoring of T and u (cf. fig. 1 for an illustration). If we are solving eq. (2.2), and notice that T → T − c (respectively T → T + c ), then we need to go over into eq. (2.7), with the initial condition u = 0 (respectively u = 1). If we are solving eq. (2.7), and notice that u → 0 (u → 1), then we need to go over into eq. (2.2), with the initial condition e r = e r (T − c ) (e r = e r (T + c )). It is possible that the system enters and exits the mixed phase from the same side (for instance, if T max = T c ), or from different sides (if the transition is passed through on the way towards higher or lower temperatures).

Non-perturbative thermodynamic functions for a radiation plasma
As an essential ingredient to eqs. (2.2) and (2.7), we need the thermodynamic energy density and pressure of the radiation plasma, e r and p r . These are often parametrized through degrees of freedom g * or h * , as e r ≡ g * π 2 T 4 /30 and e r + p r = T s r ≡ h * 2π 2 T 4 /45, where s r is the entropy density. If the plasma is very weakly coupled, g * is to a good approximation constant and h * ≃ g * , however our focus is on self-interacting plasmas, where the interactions can become strong as well. 3 In the latter case, g * and h * decrease rapidly at low temperatures, and their complete functional forms are needed.
It turns out to be convenient to parametrize the thermodynamic information through the entropy density, s r . On one hand, this is because we need s r for eq. (3.2); on the other, because s r can be precisely studied with lattice simulations (cf. refs. [9,10] and references therein). Denoting by T c the critical temperature, and setting N c = 3 from now on, the results of the deconfined phase of a Yang-Mills plasma can be represented as [10] The determination of s r /T 3 is more difficult in the confined phase, as the results soon become exponentially small. Fitting to the tabulated results from ref. [10], 5 viz. which appear to be consistent with ref. [11], we model the low-T region with the ansatz In this section and in sec. 2.5, we make use of non-perturbative information, so that the coupling can be arbitrarily strong, whereas for the friction coefficient discussed in sec. 2.4, reliable non-perturbative information is not available. Then we extrapolate weak-coupling predictions, strictly speaking only applicable for α < 0.3, to a strongly coupled regime, introducing an estimate of the corresponding error along the way. 4 Ref. [10] gives T = 3.433T c as the transition point between the two functional forms; we have replaced this with the value at which the curves cross each other, with the approximate coefficients at our disposal. 5 The first number can be found in the text, not the table. The last number appears to contain a typo in the table of ref. [10]; we have reconstructed the correct value from the e r and p r given on the same line.
The transition is of the first order, so that s r /T 3 displays a discontinuity at T = T c . For the conversion between T c and Λ IR , we estimate T c ≃ 1.24Λ IR [12]. Given s r = dp r /dT , the other thermodynamic functions can be obtained as Furthermore, in order to evaluateė r =Ṫ c r , we need the heat capacity c r = ∂ T e r = T ∂ T s r .
To keep c r continuous, we have moved the matching point in eq. (2.8) to T = 4.863T c when evaluating ∂ T s r . At low temperatures, in turn, eq. (2.11) implies Restricting to the SU(3) plasma is a special case, however this is the system for which the most reliable non-perturbative information is available. In addition, it entails a weak first-order transition, which is typical of many other thermal systems.

Perturbative thermodynamic functions for a thermalized inflaton
As the system heats up, the inflaton field might equilibrate as well (see, however, the discussions in secs. 2.7 and 3.2). Around the minimum of the potential, the inflation is a weakly coupled massive scalar field, whose interactions are suppressed by powers of 1/f a . Then the effective potential V contains a temperature dependent part, which contributes to the thermodynamic functions p ϕ and e ϕ according to eq. (2.3).
The starting point for the evaluation of p ϕ and e ϕ is the 1-loop expression for a thermal effective potential, We omit the T -independent vacuum part of V eff in the following. Changing variables to a form convenient for a numerical evaluation, the pressure and energy density from eq. (2.3) then obtain the contributions where n B (ǫ) ≡ 1/(e ǫ/T −1) is the Bose distribution. The tree-level potential V 0 should contain no T -dependence, so that we have set V 0,T to vanish, however this requires some discussion of thermal mass corrections, to which we return in sec. 2.5. Finally, eqs. (2.2) and (2.4) contain the contribution of ϕ to the heat capacity, In the massless limit, i.e. m ≪ πT , eqs. (2.16)-(2.18) amount to the substitutions g * → g * + 1 and h * → h * + 1 in the number of effective degrees of freedom. 6

Friction coefficient
A key role for the heating-up dynamics, according to eq. (2.2), is played by the friction coefficient Υ. If Υ = 0, like in standard cold inflation computations, there is no source term for the temperature evolution, and any possible initial temperature just redshifts away. It has been realized, however, that the assumption Υ = 0 is mathematically troublesome. The problem is that even if the T = 0 solution represents a fixed point, it can be an unstable one. Just a small perturbation may drive the system to another fixed point, where T > 0 and Υ > 0 [13,14]. The properties of this thermal fixed point constitute the topic of sec. 3.1.
In general, Υ is a function of the frequency, ω, at which the system is probed [15]. Then the full equation of motion does not have a local form, but rather contains a dispersive integral over the medium response.
A local evolution equation is obtained around the global minimum, where the frequency scale can be replaced by the corresponding mass scale, ω → m. Before the system settles to the global minimum, the situation may be intuitively probed by replacing the frequency scale by the curvature of the potential, ω → ω dyn ≡ max(0, V ϕϕ ) [16]. Then, if V ϕϕ ≤ 0, temperature is the only scale at early stages of inflation. But, as mentioned, T = 0 is an unstable fixed point in this setup. As the focus of the current study is the heating-up period, we will adopt the replacement ω → m throughout, with the understanding that at early stages of inflation this is just a recipe.
Like the thermodynamic functions in sec. 2.2, the determination of Υ requires lattice simulations. Lattice simulations come in two different variants. At very high temperatures, T ≫ Λ IR , so-called classical real-time simulations can be employed, and this is the method used for estimating Υ IR in eq. (2.20) [17], as well as the shape in eq. (2.19) [18]. However, when T < ∼ Λ IR , these effective-theory type setups should be replaced by full four-dimensional lattice simulations. Unfortunately, extracting real-time information from the latter is exponentially hard (cf., e.g., refs. [19,20] been launched [21,22]. For this reason, our estimates contain a systematic error, reflected by the x-dependence introduced in eq. (2.21). Now, for a qualitative understanding, it is often sufficient to consider limiting "thermal" (ω ≪ T ) or "vacuum" (ω ≫ T ) frequencies [15]. 7 Then where d A ≡ N 2 c − 1 and κ ≃ 1.5. The latter represents the vacuum decay width for the process ϕ → gg. For the gauge coupling, we adopt a leading-logarithmic running value, representative of a Yang-Mills plasma, The first term in the square root serves as an (arbitrary) infrared (IR) regulator, so that any value of T or ω can be inserted; we will check the IR sensitivity of the results by varying the parameter x in the range x ∈ (0.2...2.0). Nevertheless the expression is guaranteed to be physically meaningful only for max{2πT, ω} ≫ Λ IR , so that α ≪ 0.3.

Mass correction
Apart from a friction coefficient, the Yang-Mills plasma in general induces a mass correction to the inflaton field [15]. This is again a function of the frequency, ω. In principle we could carry out a discussion similar to Υ, adapting a weak-coupling computation from T ≫ Λ IR [16] to a strongly coupled regime through a modelling of α. However, partial non-perturbative lattice information is also available, so we would like to make use of it, even if it is not exactly what is needed. Let us explain the issues. The lattice simulations are usually viewed in the context of the QCD axion mass. As a rule, it is (implicitly) assumed that the axion has no mass at tree level, so that all of it is generated by the SU(3) gauge dynamics. Then, the mass correction can be evaluated at ω = 0, in which case it is proportional to the so-called (Euclidean) topological susceptibility. Though the problem is technically challenging, results have become available (cf., e.g., refs. [23,24] and references therein), and we return to them presently (cf. eq. (2.22)).
However, a mass determined at ω = 0 is correct only if there is no "bare mass" from an UV theory. This is a problem, since our inflaton potential V 0 already contains a mass; it 7 In the numerical solutions, we use the full interpolating function as estimated in ref. [18], viz.
is envisaged to have been generated by an UV gauge theory, at a scale higher than the IR one on which we focus. In this situation, the contribution to the mass through the SU (3) topological susceptibility is only a correction, and its determination at ω = 0 rather than ω ≃ m represents an uncontrolled approximation from the physics point of view. Despite these reservations, let us estimate how large the mass correction could be. We denote by χ topo the SU(3) topological susceptibility, and by t 0 an auxiliary quantity often used for setting the scale in lattice simulations. Then, at T = 0, t 2 0 χ topo = 6.67(7)×10 −4 [23], whereas examples of thermal values are t 2 0 χ topo = 2.25(12) × 10 −5 at T 8t 0 = 1.081 and t 2 0 χ topo = 3.43(27) × 10 −6 at T 8t 0 = 1.434 [24]. Inserting a conversion to the critical temperature, T c √ t 0 = 0.2489(14) [12], a rough qualitative represention, incorporating an expected functional dependence at higher temperatures, is and the corresponding mass correction from the IR gauge theory evaluates to (2.23) Let us connect δm 2 IR to the parameters that appear in V 0 . We have referred to V 0 as the tree-level potential, but let us now assume that it also includes those radiative and thermal corrections which do not change the shape of V 0 (in contrast, shape-changing structures lead to what we have denoted by V , as discussed in sec. 2.3). Let then m 2 0 be the value of the mass before the inclusion of the IR contribution. Then, at zero temperature, whereas a would-be thermal mass squared reads According to eq. (2.22), δm 2 IR,T < δm 2 IR | T =0 , so the thermal correction in eq. (2.25) is negative: the effective mass squared decreases at high temperatures.
After these order-of-magnitude estimates, let us explain why the mass correction should be unimportant. First of all, inserting numerical values from eq. (2.27) and estimating the bare mass as m 2 ∼ Λ 4 UV /f 2 a , in accordance with eq. (2.23), the scale of the UV gauge theory is Λ UV ∼ mf a ∼ 10 −3 m pl . Now, according to sec. 3.2, the solutions fall in two different classes. If T max ≥ T c , then |δm 2 i.e. the thermal mass correction is exceedingly small. If T max < T c , the system stays in the confined phase, and there is no thermal mass correction.
It is for these reasons, as well as the conceptual issues explained at the beginning of this section, that we omit thermal mass corrections in the following.

Inflaton potential and parameter choices
In order to study heating-up in a semi-realistic framework, we consider axion-like (or natural) inflation [4], and fix the parameters of the potential to agree with Planck data [5]. For this purpose, we adopt the predictions of cold inflation, despite the fact that the solution is technically unstable (cf. sec. 2.4). The coupling α does not appear in these predictions, and can thus be freely varied for understanding heating-up dynamics. In principle it would be interesting to verify a posteriori how significantly phenomenological predictions are altered by thermal effects if high temperatures are reached already during the inflationary stage (for a review of warm inflation see, e.g., ref. [25]). However, this dynamics is physically distinct from and takes place much earlier than the heating-up stage that we are interested in.
For the inflaton potential, we take the ansatz whereas V 0,T is approximated as small (cf. sec. 2.5). In the numerical estimates the parameters are fixed to benchmark values from ref. [16], where t ref denotes the time at which we start the simulation, conveniently chosen as .
(2.28) This is within O(1) of the initial inverse Hubble rate.

When is the thermalization assumption self-consistent?
To conclude this section, let us return to the important question of when the temperature is a useful concept. An upper bound extending almost up to the Planck scale was presented in sec. 1, however it was based on the assumption that the Hubble rate is already dominated by the energy density carried by radiation. The question of how fast a general system equilibrates is a hard one. Frequently, the dynamics following inflation is studied by solving classical field equations of motion. 8 The problem is that classical field dynamics can be correct only for large occupation numbers, not in the typical domain where the occupation is of order unity. But it is precisely momenta from the latter domain, p ∼ πT , which carry most of the radiation energy density. In other words, the issue of thermalization cannot be properly resolved with classical field theory. In the heavy-ion context, where thermalization of non-Abelian systems has been studied extensively, the method of choice relies nowadays rather on effective kinetic theory [28]. If we do think in the language of effective kinetic theory, we can draw diagrams responsible for thermalization. Examples for a non-Abelian plasma and for the inflaton, respectively, are illustrated in fig. 2. The gauge plasma equilibration rate (i.e. the thermally averaged amplitude squared) is then Γ g ∼ α 2 T within the weak-coupling expansion, whereas the inflaton equivalent is Γ ϕ ∼ Υ IR . In the subsequent sections, we will compare these with the Hubble rate a posteriori, while the computations themselves are carried out in the presence of a temperature-like parameter, reminiscently of the scenario of warm inflation [25]. Here it is appropriate to remark that the general objections to warm inflation, raised in ref. [29], are avoided in the non-Abelian axion inflation context, as it is possible to have a large thermal friction coefficient (cf. sec. 2.4) without inducing a large thermal mass correction (cf. sec. 2.5).

Temperature evolution
Before attacking the full numerical solution of eqs. (2.1), (2.2), and (2.7) (cf. sec. 3.2), we introduce the concept of a stationary temperature, T stat , which already helps us to understand the parametric dependence of the temperature scale reached (cf. sec. 3.1).

Stationary temperature as a qualitative estimate
The existence of a stationary temperature at intermediate stages of warm natural inflation follows from the argumentation in refs. [13,14]. Physically, this corresponds to a situation in which the energy released from the inflaton to radiation through friction, precisely balances against the energy diluted by the Hubble expansion. After a while, T starts to increase above T stat , obtaining a maximal value, T max . While estimating T max requires the solution of a coupled set of differential equations, it is easier to determine T stat , as the equations are algebraic, and we may furthermore employ slow-roll approximations in them. Yet, as consolidated by the numerical studies in sec. 3.2, T stat already gives an order-of-magnitude estimate of T max .
In eq. (2.2), we search for a stationary solution, withė r − TV T ≃ 0. Recalling the thermodynamic relation e + p = T s, where s is the entropy density, yields the master relation As further simplifications, the Hubble rate can be approximated as H ≃ 8πV /(3m 2 pl ) during the slow-roll period, and we may furthermore setφ →φ(t ref ) in V and V ϕ .
We illustrate the solution originating from eq. (3.2) in three qualitatively different cases: (i) confining plasma, non-thermal inflaton. In this case we insert the equation of state from sec. 2.2, s → s r , and assume that the inflaton does not thermalize. This assumption can often be justified a posteriori, cf. sec. 3.2.
(ii) confining plasma, thermalized inflaton. We add the contribution of a thermalized inflaton to the entropy density, s → s r − V T , as specified in sec. 2.3.
(iii) confining plasma with a light degree of freedom and thermalized inflaton. We further add one free massless bosonic degree of freedom to the radiation plasma, in order to illustrate the qualitative influence of a dark photon or a dark light pion. Then s → s r − V T + 2π 2 T 3 stat /45 in eq. (3.2).
For the three cases defined, the left and right-hand sides (lhs, rhs) of eq. (3.2) are illustrated in fig. 3(left). The resulting values of T stat /m pl , from the crossings of the respective curves, are plotted in fig. 3(right).
It can be observed from fig. 3(right) that there is a specific domain of Λ IR at which the behaviour of the system changes. Let us denote the smallest value of Λ IR for which the system stays stationary at exactly the critical temperature by [ Λ 0 ] stat . This corresponds to A similar logic will be applied in sec. 3.2, with however the stationary temperature replaced by the maximal one, yielding then an outcome denoted by Λ 0 . In between [ Λ 0 ] stat and Λ 0 , the system heats up to above T c , and experiences two nearby phase transitions, before and after this moment. At Λ IR > Λ 0 , T max < T c , and no phase transition takes place. Now, in terms of fig. 3(left), we find that when Λ IR = [ Λ 0 ] stat , then the lines cross in a domain where the rhs curve is flat. The cusp in the rhs curve is where the behaviour of Υ changes, from Υ UV at low temperatures to Υ IR at high temperatures. Therefore, [Λ 0 ] stat is determined by Υ UV , and independent of the non-perturbative physics of Υ IR that originates from sphaleron dynamics. The latter is important in the regime Λ IR ≫ [Λ 0 ] stat .

Maximal temperature from a numerical solution
We now move on to consolidate the qualitative results from sec. 3.1 through a full solution of eqs. (2.1), (2.2), and (2.7). For this we consider the same three cases as defined in sec. 3.1 and illustrated in fig. 3. As initial conditions, we take T (t ref ) ≃ 0.2T stat ,φ(t ref ) from eq. (2.27), and the slow-roll evolution rateφ(t ref ) ≃ −V ϕ /(3H). However, the solution is an attractor, and therefore soon independent of the initial conditions. The domain in which the system heats up exactly to T c is denoted by Λ 0 , with the numerical value Λ 0 ∼ 2 × 10 −8 m pl .
Examples of solutions are shown in fig. 4, and a scan of the resulting values of T max in fig. 5. We conclude that • if Λ IR < Λ 0 , T max > T c , and the system undergoes a phase transition as it cools down.
For [Λ 0 ] stat < Λ IR < Λ 0 , it also undergoes a nearby previous transition as it heats up, as visible in fig. 4(left). However, in this domain α 2 T ≪ T ≪ H, so it is questionable whether the temperature has a literal meaning during the first transition (cf. sec. 2.7).
• if Λ IR > Λ 0 , T max < T c , and the system undergoes no phase transition. Nevertheless it heats up to a high temperature, T max ∼ Λ IR , if the plasma is confining. If the inflaton thermalizes, the maximal temperature is cut off by the inflaton mass, T max ∼ min(m, Λ IR ). However, it appears unlikely that the inflaton thermalizes, because its would-be thermalization rate ∼ Υ IR is much below the Hubble rate, cf. fig. 4(right). If the plasma includes g * massless degrees of freedom, the maximal temperature stays at T max ∼ Λ 0 (2d A /g * ) 1/4 , irrespective of the value of Λ IR . However, this consideration assumes that the massless degrees of freedom thermalize, i.e. that their full phase space can be filled to accommodate the entropy released from inflaton oscillations. 9 It can be observed from fig. 5 that in cases (ii) and (iii), there is a large uncertainty in the value of T max . This originates from the gauge coupling α, through the parameter x (cf. eq. (2.21)). The reason is that these would-be solutions lie deep in the confined phase, where our treatment of Υ is not reliable. However, as discussed above, the physical significance of these solutions is questionable for another reason as well, namely that the thermalization assumption is hard to consolidate. It is a lucky coincidence that both the technical and conceptual uncertainties can be reflected by the same error bands. Our setup is self-consistent only for Λ IR < Λ UV ∼ 10 −3 m pl , so we restrict the axis to this domain. The grey bands show the effect of varying the parameter x in eq. (2.21) in the range x ∈ (0.2, 2.0). The uncertainties are huge for cases (ii) and (iii) at large Λ IR , as we then need to evaluate the friction coefficient Υ deep in the confined phase. By T c we denote the critical temperature.

Physical implications for gravitational waves
The patterns observed in sec. 3.2 have a number of potential implications for primordial gravitational waves, to which we now turn.

Total energy density in the gravitational wave background
A thermal plasma necessarily generates a spectrum of gravitational waves, whose energy density contributes to the effective number of massless degrees of freedom, N eff [30]. As long as the maximal temperature is below ∼ 10 −2 m pl , the contribution to N eff is small [31][32][33][34][35]. At the same time, our framework is consistent as long as Λ IR < Λ UV ∼ 10 −3 m pl . Given that T max < Λ IR for large values of Λ IR , we then also have T max < Λ UV . Therefore none of the heating-up scenaria that we have found is excluded by constraints from N eff .
to the dark sector one, whereby Standard Model values can be adopted for the entropy densities [37], gives f (i) 0 < ∼ 5 × 10 3...5 Hz and f (iii) 0 < ∼ 8 × 10 6...8 Hz. A part of this range can be probed by interferometers, like again the Einstein telescope and ultimately perhaps DECIGO.
In figs. 1(left,middle), we also see another phase transition, passed as the system heats up. Concerning its significance, two issues should be raised. The first is that as the gravitational energy density scales as ∼ 1/a 4 , the signal from the first transition is diluted by a factor ∼ e −4∆N compared with the second one, where ∆N is the number of e-folds between the transitions. If ∆N ≫ 1, then the signal gets diluted away. Only for a fine-tuned value Λ IR ∼ 1.5 × 10 −8 m pl , when the transitions are immediately adjacent to each other, could the dilution be less spectacular. Second, as visible in figs. 4(left,middle), the Hubble rate far exceeds the temperature during the first transition, and also during the second one if the two transitions are nearby. As the thermalization rate is ∼ α 2 T ≪ H (cf. sec. 2.7), it is questionable whether the temperature is a physically meaningful notion in this case. Despite these concerns, it could be that non-equilibrium fluctuations produced during the first epoch could have a physical effect, for instance by serving as nucleation seeds which would permit for the second transition to proceed in a non-standard manner (cf., e.g., ref. [38]).

Summary and outlook
The purpose of this study has been to estimate the maximal temperature that strongly coupled dark sectors may reach. The basic point is that in the confined phase of such theories, thermodynamic functions such as the entropy density and the heat capacity are exponentially small. Therefore, even a small release of energy density from the inflaton field can heat up the system by a large amount. For the very largest confinement scales, Λ IR ∼ (10 −8 − 10 −3 )m pl , we find that the system heats up to close to the critical point, even if it remains just slightly below T c (cf. fig. 5). Interestingly, this most interesting scenario is treated most reliably by our methods, as the gauge field thermalization rate clearly exceeds the Hubble rate, so that there is no doubt about the validity of temperature as a physical notion (cf. fig. 4(right)).
As a particular consequence of such dynamics, we have considered the gravitational wave background produced by thermal fluctuations around the heating-up epoch. In the frequency window (10 −4 − 10 2 ) Hz, relevant for the LISA, Einstein telescope, and DECIGO interferometers, we predict a background increasing monotonically as ∼ f 3 0 , with a coefficient proportional to the maximal shear viscosity of the plasma phase (cf. eq. (4.1)). At the highest frequencies, the signal could be marginally observable in the future, though quantitative sensitivity studies would be needed for confirming or refuting this prospect.
A complementary consequence is reached if the confinement scale is lowered, Λ IR < 10 −8 m pl . Then the system heats up above the critical temperature, confirming the possibility of a firstorder phase transition in a dark sector [36].
We should underline that we have on purpose kept our study on a rather general level, with the hope that it is then also more broadly applicable. Adding specific model assumptions, further issues could be addressed. Notably, the value of the maximal temperature, and in particular whether it is above T c or not, has implications for dark matter production, but the details are very model-dependent (cf., e.g., refs. [1][2][3]).
Apart from phenomenological issues, there are also theoretical ingredients that could be built into our framework. Hoping that exploratory low-temperature investigations of the friction coefficient [21,22] turn ultimately into a semi-quantitative tool, the error bands in fig. 5 could be reduced. Were the same to happen with the shear viscosity η, the gravitational wave estimate in sec. 4.2 could be sharpened. Inserting assumptions about the couplings of the inflaton to both the dark and the visible sector, and of the sectors between each other, two different temperatures could be tracked. Finally, investigating the non-equilibrium physics of a system possessing two adjacent transitions (cf. fig. 1(middle)) might reveal interesting gravitational wave signatures. In the last case, it should be recalled that the transition takes place during a period in which inflaton oscillations dominate the energy density, resulting in changes to the normal predictions that assume a radiation-dominated universe.