Gravitational waves from scale-invariant vector dark matter model: probing below the neutrino-floor

We study the gravitational waves (GWs) spectrum produced during the electroweak phase transition in a scale-invariant extension of the Standard Model (SM), enlarged by a dark U(1)D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ U(1)_{D} $$\end{document} gauge symmetry. This symmetry incorporates a vector dark matter (DM) candidate and a scalar field (scalon). Because of scale invariance, the model has only two independent parameters and for the parameter space constrained by DM relic density, strongly first-order electroweak phase transition can take place. In this model, for a narrow part of the parameter space, DM-nucleon cross section is below the neutrino-floor limit, and therefore, it cannot be probed by the future direct detection experiments. However, for a benchmark point from this narrow region, we show the amplitude and frequency of phase transition GW spectrum fall within the observational window of space-based GW detectors such as eLISA.


Introduction
The detection of GWs [1] has opened up a new and independent avenue for probing of dark matter [2]. These waves are ripples in the fabric of space-time generated by energetic and violent sources such as black hole and neutron star binaries, extreme mass ratio inspirals, and first order cosmological phase transitions. The main targets of ground-based GW detectors are black hole and neutron star binaries with best sensitivity at frequencies O(10 2 ) Hertz, while space-based detectors are most sensitive to milli-Hertz or deci-Hertz frequencies [3]. The stochastic background of primordial GWs produced during first order electroweak phase transition is a physical sources of GWs in this frequency band [4].
Cosmic phase transitions occur when the temperature drops below a critical temperature leading to the transition of the Universe from a symmetric phase to a phase of broa e-mail: mohamadnejad.a@lu.ac.ir (corresponding author) ken symmetry (for a recent review see [5]). In the SM, electroweak phase transition, as well as QCD phase transition, is of second order [6,7] and does not generate the GW signal. However, early Universe might be in a state where not only the gauge symmetry was present but also a (classical) scaleinvariant or conformal symmetry was realized, preventing any massive parameters in the Lagrangian. Then the quantum effects must have broken the scale-invariant symmetry, generating a nonzero VEV of the scalar field(s) and all the mass terms of massive particles via Coleman-Weinberg mechanism [8]. Particularly, conformal DM models with Higgs portal are attractive because they can solve DM problem and at the same time can generate a strongly first order phase transition [9][10][11][12]. GWs due to first order phase transition have been studied within models where the scale-invariant symmetry is broken due to Coleman-Weinberg mechanism [13][14][15][16][17][18][19][20][21] or models with DM candidate [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41]. Conformal symmetry also proposed as a possible solution for hierarchy problem [42].
Strongly first order electroweak phase transition can take place in the early Universe when two local minima of free energy (potential) co-exist for some range of temperatures. It is a necessary condition in generating the observed baryon asymmetry in the universe and provides one of the three Sakharov conditions [43], i.e., an out-of-equilibrium environment, in the framework of electroweak baryogenesis (for a recent review see [44]). On the other hand, this violent phenomenon can lead to large anisotropic fluctuations in the energy-momentum tensor generating stochastic GW background [45][46][47][48][49][50]. This signal can potentially be probed in future space based GW detectors such as Laser Interferometer Space Antenna (LISA) [51,52], Big Bang Observer (BBO) [53], Deci-hertz Interferometer Gravitational wave Observatory (DECIGO) [54], and Ultimate-DECIGO (U-DECIGO) [55].
After two local minima of the free energy co-exist at a critical temperature, the relevant scalar field can quantum mechanically tunnel into the new phase. This phenomenon continues via the nucleation of bubbles which expand and eventually collide with each other leaving a significant background of GWs. According to simulations of GW backgrounds from cosmic phase transitions, there are three spectral contributions: (1) the collision spectrum which is the direct result of bubbles of true vacuum colliding [56][57][58][59][60][61], (2) sound wave spectrum which is the consequence of the fluid dynamics following such collisions, dominating in most relevant scenarios [62][63][64][65], and (3) the turbulence spectrum, which is generally subdominant [66][67][68][69][70][71]. GWs can prevail to the present times and perhaps be detected in space-based GW detectors.
Scale-invariant extensions of the SM can also provide a DM candidate (for recent papers see e.g. [72][73][74][75][76][77][78][79][80][81][82][83][84]). Particle nature of the DM is another important puzzle in particle cosmology [85]. There are in principle three ways to search for such exotic particle: (1) direct detection, (2) indirect detection, and (3) collider searches. Direct detection experiment such as the LUX [86], PandaX-II [87] and XENON1T [88] are gradually approaching the neutrino backgrounds which is usually considered as the ultimate sensitivity of future direct detection experiments [89]. Neutrino background (neutrinofloor) puts a limit on discovery potential of DM. To this day, DM search experiments have not found any evidence. However, null results of DM detection does not exclude the possibility of observing a GW signal from a dark sector. And in some models GW signals could be a unique probe of the thermal DM paradigm. In the absence of DM signal below neutrino-floor, GW experiments may serve as a new approach to probe the DM models. Consequently, GW detectors can be a vital tools in exploring possibilities for DM models, complementing existing efforts at colliders, direct and indirect detection experiments.
In this paper, we will study a conformal model [12] which is well motivated from a DM perspective or by naturalness arguments. The model is a scale-invariant extension of the SM, enlarged by a dark U (1) D gauge symmetry which provides a viable vector DM candidate. There are only two additional fields as well as free parameters in this model and the model can overcome constraints such as DM relic density, and direct and indirect detection upper bound limits. For the parameter space constrained by DM relic density, strongly first-order electroweak phase transition can also take place. The main interest of the present paper is the GW signal produced during the electroweak phase transition. We use a benchmark point of the parameter space with below the neutrino-floor DM-nucleon cross section and show that the amplitude and frequency of phase transition GW spectrum fall within the observational window of eLISA. Since, this particular choice of the parameter space is not constrained by colliders, and direct or indirect detection, therefore, GW signal plays an important role in probing the model for the chosen benchmark point. This paper is structured as follows. In Sect. 2 we introduce the model and review DM phenomenology. In this section, we choose a benchmark point below the neutrino-floor for the rest of the paper. In Sect. 3, we study effective potential. GW signal produced during first order electroweak phase transition as well as its discovery prospects are presented in Sect. 4, after which Sect. 5 comprises a summary and our conclusion.

Review of the model and DM phenomenology
Let us first give an overview of the model presented in [12]. The beyond SM fields content of the model are a complex scalar field (φ) and a vector field (V μ ). The model is a conformal extension of SM, enlarged by a dark Abelian gauge symmetry. The scalar field φ has a unit charge under dark U (1) D symmetry and V μ is the corresponding Abelian gauge field serving as DM particle. These two fields are neutral under SM gauge group and SM fields are singlet under U (1) D . There is also a discrete Z 2 symmetry, under which SM particles are singlet and the vector field V μ and the scalar field φ transform as follows: Note that, Z 2 symmetry forbids the kinetic mixing between the vector field V μ and SM U Y (1) gauge boson B μ , therefore, the vector field V μ is stable and can be considered as a viable DM candidate. The conformal model may be organized into three sectors: (1) the visible sector which consists of the SM fields without Higgs potential (L V S ), (2) the dark sector which consists of the vector DM, V μ , together with the scalar field φ (L DS ), and (3) scale invariant tree-level potential (V tree ). The Lagrangian is given by where The most general scale-invariant potential which is renormalizable and invariant under Z 2 and gauge symmetry is where H is the Higgs doublet. In Eq. (4), the third term is the only connection between the dark and the visible sector.
In unitary gauge, H † = 1 √ 2 (0 h 1 ) and φ = 1 √ 2 h 2 , therefore, tree-level potential is given by where h 1,2 are real scalar fields. Vacuum stability requires λ H,φ > 0 and λ φ H < 0,. Furthermore, non-zero VEV of h 1,2 scalar fields demands The Local minimum of the two-variable potential (5) defines a direction in field-space known as flat direction [90]. Along this direction V tree = 0, while in other directions V tree > 0. Therefore, tree-level potential only vanishes along the flat direction where Naturally, we expect higher-loop contributions be dominated along this direction and they should determine the local minimum of the full potential containing higher-loop effects. Indeed, for some mass spectrum of the model, 1- , gives a small curvature in the flat direction with V 1−loop e f f < 0. Therefore, considering 1-loop effect, the potential has a global minimum point in field space and consequently, symmetry breaking can take place. We assume ν 1 and ν 2 are VEVs of h 1 and h 2 where ν 1 = 246 GeV. Now consider mass eigenstates h and ϕ, where α is the angle between flat direction and h 2 axis in field-space. Therefore, along the flat direction h = 0, and all massive particles get mass when ϕ = 0. In our formulation, h is perpendicular to the flat direction and we identify it as the SM-like Higgs observed at the LHC with m h = 125 GeV. At the classical tree-level, the scalon field ϕ is massless, however, the 1-loop corrections give a mass to this field via Gildener-Weinberg mechanism [90]. Regarding 1-loop effect, the scalon mass is given by where m V,W,Z ,t being the masses for vector DM, W and Z gauge bosons, and top quark, respectively, and ν = ν 2 1 + ν 2 2 . After symmetry breaking, all the four dimension-less parameters of the model, i.e., λ H,φ H,φ and g v , will be determined by DM mass m V , m h , ν 1 , and ν 2 . Since we have already determined m h and ν 1 , therefore, the model has only two independent parameters. Here, we choose, DM mass m V and ν as our two-dimensional parameter space. According to Eq. (8), m ϕ is also determined by this two-dimensional parameter space. Fig. 1 DM-nucleon cross section for a parameter space already constrained by DM relic density, h 2 = 0.12 [91] We have recently studied the DM phenomenology of such a model [12]. Here we briefly report the main results. In our model, for a narrow region of the parameter space, a direct detection of the vector DM is hopeless due to the neutrino-floor which represents an irreducible background, see Fig. 1.
In this figure, spin-independent (SI) DM-nucleon cross section versus DM mass is depicted. PandaX-II [87] upper bound as well as neutrino-floor [89] limit is shown for comparison. As it is seen, for a narrow region of parameter space, DM-nucleon cross section has a dip below the neutrinofloor. The origin of this dip, is the proportionality of DMnucleon cross section with ( 1 [12]. Therefore, for the parameter space around m ϕ m h , we expect such a dip.
As we mentioned before, only a small portion of the parameter space is below neutrino-floor. This region of parameter space is beyond the sensitivity of future Direct or indirect detection experiments [12]. However, we show that one can probe this blind-spot of the model using space-based GW detectors. Since the region below the neutrino-floor is narrow, we choose only one benchmark for the rest of this paper, see Table 1.

Effective potential
The effective potential is composed of three parts: (1) classical or tree-level potential, (2) zero temperature 1-loop potential known as Coleman-Weinberg potential, and (3) 1-loop finite temperature potential. In the following we study these three pieces of full effective potential. Table 1 A benchmark point of the model and corresponding DM relic density, SI DM-nucleon cross section, and DM total annihilation cross section. For this benchmark point, DM-nucleon cross section is below the neutrino-floor limit In the white region, M 2 ϕ will be complex

Tree-level potential
The tree-level potential is given in Eq. (5). In the scalar sector, tree-level field-dependent masses correspond to the eigenvalues of the Hessian matrix: and are given by For the top quark and the gauge bosons the field dependent masses are where λ t denotes the top quark Yukawa coupling, and g v , g, and g are dark U (1) D , SU (2) L and U (1) Y gauge couplings, respectively. According to Eq. (11), it is obvious that M 2 t , M 2 V , M 2 W , and M 2 Z are positive. As Fig. 2 shows, M 2 h is also positive, however, M 2 ϕ < 0 except along the flat direction where M 2 ϕ = 0. Here, we exclude the field space with M 2 ϕ < 0 and only consider flat direction.

Coleman-Weinberg potential
The Coleman-Weinberg potential is a sum of 1PI 1-loop diagrams with arbitrary numbers of external fields and particles running in the loop and it is given by [8] where C k = 3/2 (5/6) for scalars/spinors (vectors), M k is the tree-level field-dependent mass of particle k given in Eqs. (10) and (11), and g k presents the number of degrees of freedom given by where s k is the spin, N k the number of colors and q k = 1 (2) for neutral (charged) particles.
In Eq. (12), to get a real Coleman-Weinberg potential, the field dependent mass squared of particles can not be negative. Thus the allowed field space should be along the flat direction. On the other hand, along this direction M ϕ = 0, and we do not consider scalon field contribution in Coleman-Weinberg potential (12). Moreover, the Goldstone bosons are massless along the flat direction and they do not contribute in the minimum of the tree-level potential. Therefore, we do not consider field dependent masses of Goldstone bosons in Eq. (12) as well.
Note that, along the flat direction, for five remained field dependent mass contributions, we can substitute M k → m k ν ϕ, where m k is the measured mass of particle k. Regarding this substitution in Eq. (12) results the well-known Gildener-Weinberg formula [90] V 1 GW = Aϕ 4 + Bϕ 4 ln where As we mentioned before, along the flat direction V tree = 0. Therefore, to find the true vacuum, one should find the minimum of the 1-loop potential (14), given by Combining Eqs. (14) and (16) we can substitute RG scale and find a simple expression for the 1-loop potential in terms of the true vacuum expectation value ν and B coefficient: From the above equation, one can find m 2 ϕ = Eq. (8). Although, the ϕ scalar obtains a radiatively generated mass at 1-loop level, Goldstone bosons remain massless to all orders in perturbation theory.

Finite temperature potential
Now we study finite-temperature 1-loop effective potential which enables us to compute scalar field vacuum expectation values, in the background of a thermal bath with temperature T . The 1-loop finite-temperature corrections are given by [92] with thermal functions vanishing as T → 0. Although these integrals cannot be expressed in terms of standard functions, their numerical evaluation is rather straightforward and one can approximate them in different limits. For example, in the high temperature limit (x 1), J B (x) and J F (x) have very useful closed forms (see appendix C in [92]), where a b = π 2 exp 3 2 − 2γ E , a f = 16 a b and γ E denotes the Euler-Mascheroni constant. The thermal functions J B (x) and J F (x) also have a useful expansion in terms of modified Bessel functions of the second kind (see Appendix A) The above sum representations are convergent as m → ∞ and in this limit J m (x) → J (x). In Fig. 3, we have depicted J (x), J high-T (x), and J m (x). As it is seen in this figure, with the log term included, high-T expansion for the thermal functions is accurate to better than 10 percent even for x ∼ (1 − 1.5) (depending on the function), but breaks down completely beyond that. However, the summation in Eq. (21) can be truncated at a few terms, for example m = 2 for J F and m = 3 for J B , and still yield a very good accuracy (see Fig. 3).
In order to consider as well the resummation of the Matsubara zero modes, it is essential to resum the thermal masses by substituting M 2 → M 2 tree + . In the standard method, is taken to be the leading contribution in temperature to the 1-loop thermal mass [93]. For scalars can be estimated by differentiating V 1 T with respect to ϕ yielding ∝ T 2 . This replacing automatically contains daisy contributions to all orders in the effective potential.
Using high-T approximation, the field dependent terms in logs cancel between V 1 T and V 1 CW . The x 2 term gives an overall contribution proportional to T 2 i . If we use only the leading-order contribution to i in temperature, T 2 i will be field-independent. Therefore, only x 3 term is left, which can be entrap by adding V ring . This means, effective potential contains four terms: where V ring is daisy term [93]: Including V ring amounts to resumming the IR-divergent terms to the Matsubara zero mode propagator. It is equivalent to substitute M 2 → M 2 tree + (T ) in the full effective potential, assuming that only the thermal mass of the zero mode is relevant, which means using high-temperature approximation.
The sum in Eq. (23) runs only over longitudinal degrees of freedom of the gauge bosons and scalars. The thermal masses of the gauge bosons are given by [93] where n f = 3 is the number of fermionic generations. For scalars we have where V high−T is derived from Eq. (18) using high temperature approximation (20) with only x 2 term, i.e., The factors of the scalar couplings in Eq. (25) are different from the results in the literature (e.g. Ref. [93]) since we did not include the impact of the Goldstone bosons (as the same as Ref. [20]). Considering Eq. (22) along the flat direction (where V tree = 0), we can obtain a one-dimensional effective potential, V e f f (ϕ, T ) which contains Gildener-Weinberg potential (17) and thermal contributions (18) and (23). In the next section using this potential, we study the phase transition and GW signal.

Electroweak phase transition and GW signal
The electroweak phase transition takes place after the temperature of the universe drops below the critical temperature (see Fig. 4) and the minimum with non-zero scalon VEV becomes the global minimum. For the parameters in Table 1, the critical temperature is T c = 339 GeV at which the effective potential has two degenerate minimums.
After the minimum with non-zero scalon VEV becomes the global minimum, thermal fluctuations eventually excite the field enough to cross the potential barrier. For conformal extensions of SM, the potential barrier disappears only for T = 0, thus the electroweak phase transition in scaleinvariant models is always of first order for any finite temperature.
The decisive quantity is the temperature at which phase transition proceeds via nucleation and consequent expansion of bubbles inside of which the field is in the broken phase of the model. The rate of bubble nucleation per unit of time and volume is given by [94] (T ) where is the three-dimensional Euclidean action for a spherical symmetric bubble. The differential equation minimizes S 3 , therefore, the field ϕ obtained by solving the Eq. (29) has the largest contribution to the bubble nucleation rate. The solution of (29) starts for r = 0 close to a specified true vacuum in field space with dϕ dr = 0 and asymptotically come near to a specified false vacuum as r → ∞. The nucleation temperature T n is defined as the temperature at which the probability of one bubble nucleation per horizon volume in Hubble time approaches unity: 4 1.
The Hubble parameter as a function of temperature is given by [95] where M pl ∼ 10 19 GeV and g * ∼ 100 are the Planck mass and the effective number of relativistic degrees of freedom in the thermal plasma, respectively. One can also estimates T n by the condition S 3 (T n )/T n 140 [96]. In order to solve Eq. (29) and find the Euclidean action (28), we have used AnyBubble package [97]. The result is depicted in Fig. 5. Using condition (30), we determine the nucleation temperature, T n = 47 GeV, which is much lower than the critical temperature. Therefore, phase transition proceeds after a large amount of supercooling and transition should be very strong, consequently we expect a fairly large GW signal produced at the temperature T * . For typical phase transitions with negligible reheating, T * is approximately equivalent to the nucleation temperature T n . Now we study stochastic GW background produced by strong first-order electroweak phase transitions. The resulting contributions come from three processes: • collisions of bubble walls and shocks in the plasma, • sound waves to the stochastic background after collision of bubbles but before expansion has dissipated the kinetic energy in the plasma, and  Table 2 Parameters of the phase transition for the benchmark point in Table 1 T c (GeV) T n (GeV) α β / H * 339 47 24 808 • turbulence forming after the bubbles have collided.
These three processes can coexist, and each one contributes to the stochastic GW background: All of these three contributions are controlled by four thermal parameters (see Table 2): • T n : the nucleation temperature, • α: the ratio of the free energy density difference between the true and false vacuum and the total energy density, where ρ * is given by • β : the inverse time duration of the phase transition, • v w : the velocity of the bubble wall which is anticipated to be close to 1 for the strong transitions [98].
The electroweak phase transition proceeds by the nucleation and expansion of bubbles of the new phase. Note that isolated spherical bubbles can not be a source of GWs and these waves arise during the collision of the bubbles. The collision contribution to the spectrum is given by [61] where S coll parametrises the spectral shape and is given by where Hz.
Bubble collision produces bulk motion in the fluid in the form of sound waves which generates GWs. This is the dominant contribution to the GW signal which is given by [65] The spectral shape of S sw is where Hz. (41) Bubble collisions can also generate turbulence in the plasma which its contribution to the GW spectrum is given by [70] Table 2). The black line denotes the total GW spectrum, the yellow line the contribution from sound waves, the blue line the contribution from bubble collisions, and the red line the contribution from turbulence. Sensitivity curves of the four configurations of eLISA detector (C1-C4) are also depicted [4]. The amplitude and peak frequency of phase transition GW spectrum fall within the observational window of eLISA where Hz. (44) In Eq. (43), h * is the value of the inverse Hubble time at GW production, redshifted to today, h * = 1.65 × 10 −5 T n 100 g * 100 In computing GW spectrum we have used [4,59] where κ, κ v , and κ turb denote the fraction of latent heat that is transformed into gradient energy of the Higgs-like field, bulk motion of the fluid, and MHD turbulence, respectively. The result is depicted in Fig. 6. As it is seen in this figure, the peak frequency of GW is about 0.005 Hz and it is strong enough to be detected by eLISA detector.

Conclusion
Dark matter is believed to be five times as predominant as visible matter and its gravitational impacts are seen throughout the Universe. However, we have yet to see GWs caused by dark matter, and we can think of numerous ways this might happen. Indeed, detecting such signals can constrain dark matter models. Currently, the connection between DM and GWs is largely unexplored, and much remains to be done to make existing constraints better and to study the prospects for identifying DM using upcoming experiments. In this paper, this connection is studied through electroweak phase transition. The GW spectrum from a first order phase transition is assumed to follow a broken power-law. However, the Standard Model does not feature a first-order phase transition, therefore, observed GW background of this kind would point uniquely to new physics. Detecting GWs of cosmological origin is not an easy task. There is an enormous foreground due to astrophysical sources which in principle makes detection challenging. Once the signals from merging neutron stars and stellar mass black holes have been recognized and removed, the main sources of foreground signals are galactic and extragalactic binaries. The galactic background generated by binary stars in the Milky Way is many times larger in amplitude than both the extragalactic foreground and eLISA's design sensitivity. Furthermore, the galactic background, being mostly concentrated in the galactic plane, can be removed duo to its anisotropy. Here we study a stochastic background of GWs, which is statistically homogeneous and isotopic. Therefore, the real problem is irreducible background coming from extragalactic binary stars which is dominated by emission from white dwarf binaries at a level of approximately h 2 = 10 −12 [99]. These signals typically appear in millihertz regime [100,101] and other signals of GW in this regime may be hidden by them. Nevertheless, there are some techniques that may detect stochastic GWs background in the presence of both instrument noise and foregrounds of GW spectrum, see, e.g., [102].
In this paper, we have studied the GW signatures associated with the strong electroweak phase transition in the early universe in a scale-invariant vector DM model. In this model, both DM mass and the Higgs potential originate from the effects of a real scalar field (scalon) whose mass term dynamically develop through the spontaneous breaking of classical scale-invariant symmetry. To study phase transition, we obtained effective potential including three terms: (1) treelevel potential, (2) Coleman-Weinberg 1-loop potential, and (3) finite temperature potential with daisy diagrams contributions. Since, we have considered the flat direction in the field space, the effective potential only contains the Coleman-Weinberg and thermal terms. We show that the nucleation temperature is much lower than the critical temperature, i.e., the phase transition proceeds after large super-cooling. To compute the GW spectrum, we have considered a benchmark point of the parameter space below the neutrino-floor and see that the amplitude and frequency of GW fall within the observational window of eLISA.
Although, our DM model cannot be confirmed only by detecting GWs, it certainly can be disapproved with not detecting such signals. Indeed, detecting GWs alone barely reveals the particle physics model behind the phase transition. Hence, GWs can be a vigorous tools in exploring possibilities for DM models, complementing already existing efforts at colliders, direct and indirect detection experiments. For the particular choice of the parameter space of this paper, these efforts cannot constrain the model at the present. However, with developing these experiments, our model can be distinguished from other models with the same GWs signals containing DM candidates or without them. For example, in this model, the scalon modifies the Higgs potential and, although its direct detection at the LHC is hopeless due to the small signal to background ratio, it could be possible at future collider [103]. Moreover, indirect collider searches, e.g., coming from the modification of the triple Higgs coupling or the Zh production at lepton colliders, provide another probe of the model. Therefore, both collider and GW signals can provide a realistic and complementary test of the model. However, in the absence of stronger colliders, the planned GW experiments can play an important role in investigating our model, and possibly other DM models with strong electroweak phase transition, especially in probing below the neutrino-floor.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Experimental Datasets derived from public resources and other data that support the findings of this study are available on request from the corresponding author.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Expansion of thermal functions
The natural logarithm has Maclaurin series Using the above formula in the thermal functions, If we make the substitutions we obtain