Atlas of solar hidden photon emission

Hidden photons, gauge bosons of a U(1) symmetry of a hidden sector, can constitute the dark matter of the universe and a smoking gun for large volume compactifications of string theory. In the sub-eV mass range, a possible discovery experiment consists on searching the copious flux of these particles emitted from the Sun in a helioscope setup \`a la Sikivie. In this paper, we compute the flux of transversely polarised HPs from the Sun, a necessary ingredient for interpreting such experiments. We provide a detailed exposition of photon-HP oscillations in inhomogenous media, with special focus on resonance oscillations, which play a leading role in many cases. The region of the Sun emitting HPs resonantly is a thin spherical shell for which we justify an averaged-emission formula and which implies a distinctive morphology of the angular distribution of HPs on Earth in many cases. Low mass HPs with energies in the visible and IR have resonances very close to the photosphere where the solar plasma is not fully ionised and requires building a detailed model of solar refraction and absorption. We present results for a broad range of HP masses (from 0-1 keV) and energies (from the IR to the X-ray range), the most complete atlas of solar HP emission to date.

HPs can be also produced in the early universe contributing to the dark radiation [11,20] or dark matter (DM) of the universe [21][22][23][24][25][26][27], in which case they can be searched in dark matter direct detection experiments [28] or through a small relic background of photons, which they produce because they are (slowly) decaying dark matter [21,22]. When the HP mass is above twice the electron mass, their fast decay into two electrons precludes them from being directly accessible DM candidates but they have become the prototypical force mediator in particle models where DM is self-interacting [29], which has triggered enormous attention in the intensity frontier of particle physics [30]. These particles appear as gauge bosons of new U(1) symmetries, which are the simplest extensions of the standard model of particle physics and appear copiously in its most ambitious extension, string theory [31][32][33][34][35][36][37][38][39][40][41][42][43][44].
In a huge region of parameter space, the most powerful laboratory to constrain the HP properties is the Sun, see exclusion regions labelled Sun-L and Sun-T in Fig. 1. Photons inside of the Sun can oscillate into HPs and leave it unimpeded, draining energy very efficiently. Since this energy must come from nuclear reactions, the temperature of the solar core tends to be hotter in models with HP emission, with the consequently higher neutrino flux from nuclear reactions and distorted sound speed profile. The agreement between solar models and the observations of the flux of Boron neutrinos can be used to constrain the emission of low mass bosons [45] and in particular the HP and its parameters [46]. Even stronger constraints arise from a recent global fit of helioseismology data and neutrino fluxes with realistic solar models perturbed by HP emission [47]. It is intriguing that constraints from other stelar systems, which are stronger for axions, are not as strong for small mass HPs [21,46,48].
the role of neutral atoms and metallic donors in the index of refraction for it influences strongly the conversion probability. Indeed, the presence of neutral hydrogen is responsible for displacing the locus of the resonance slightly inwards to the Sun with respect to the case in which all the plasma would be considered to be ionised. First studies along this direction were taken already some time ago and first results were presented in [69,70], which turn out to be quite accurate as order of magnitude estimates. In this paper, we ratify these estimates with a more detailed model of refraction in the Sun and exclude a possible strong influence of metals and excited states of hydrogen (with some caveats). Another complication is that, close to the photosphere, photons can travel a long distance before being absorbed, so long that the characteristics of the plasma change. In previous works, the photon-HP conversions where computed assuming that the solar plasma is homogenous in a mean-free-path but this is not the case close to the photosphere. Therefore, we also have to deal with the problem of photon-HP resonant oscillations in a non-homogenous plasma, which so far has not been discussed in this context. Perhaps, the most remarkable finding of this paper is that the total flux of HPs averaged over the region where resonant conversions happen is independent of this fact (as long as the resonance region is thin) and the previously used formulas for the Sun averaged HP emission are still correct (!). One last complication with the resonance emission is that near the surface, the Sun is not spherically symmetric. Inhomogeneities in the temperature and density of the plasma arise because the energy transport is convective and outgoing hot cells and incoming cool flows have differences in their index of refraction. The resonance region is no longer a perfect spherical cell, it becomes corrugated and timedependent. Nevertheless, the techniques developed for a spherically symmetric Sun with an effective 1D atmosphere can be applied also in this case. We find that the volume of the resonance region is slightly increased due to corrugation and the spectrum somewhat distorted because the 3D temperature profiles are steeper close to the surface and shallower inside. The result is a moderate O(1) increase of the HP flux coming from inner resonances and a decrease of the outer ones. But the resonance emission does not explain the HP spectrum in the whole energy range. Low mass HPs have resonance regions close to the solar surface where the temperature is O(eV) and emission of UV or X-ray energy HPs is exponentially suppressed. In these cases, the HP flux emitted from the bulk of the Sun dominates. At these energies, bound-bound and bound-free processes of helium and metals do contribute notably to the photon opacity and have to be included in the calculation. In this paper we use monochromatic opacities from the Opacity Project to compute this flux and find a significant contribution.
The plan of the paper is as follows. In section 2 we recall the physics of photon-HP oscillations in vacuum, homogenous and in homogenous plasmas and we derive the formulas useful for the resonant and non-resonant emission of HPs in the Sun. Section 3 is devoted to compute the ingredients needed to build a model for the refraction and absorption of photons in the solar plasma. In section 4 we present our results for the flux of transversely polarised HPs emitted from the Sun from the near infrared to the keV range, with an emphasis on the visible part of the spectrum. We first train the intuition of the reader treating the Sun as an spherically symmetric plasma and only later we discuss the implications of the convectiondriven inhomogeneities close to the surface. Finally, in section 5 we use our results to set new constraints on the HP properties and speculate on future experimental searches for the transversely polarised HP flux. Note that henceforth we will refer exclusively to transversely polarised HPs merely as HPs for the sake of simplicity. When we speak about longitudinally polarised HPs we will do so explicitly.

Hidden photons
Consider hidden photons as gauge bosons of a hidden U(1) gauge symmetry, i.e. a symmetry under which standard model particles do not transform. New particles or string excitations charged under the new U(1) and under the hypercharge of the standard model will generate kinetic mixing between the hidden photon and the hypercharge boson [2][3][4][31][32][33][34][35][36][37][38][39][40][41][42][43][44] or, at energies below electroweak symmetry breaking, with the photon and Z boson. Typically, integrating out these mediator new physics, also introduces an extra infinite tower of new operators, which in contrast to the kinetic mixing, are generically suppressed by the energy scales of the new physics (usually particle masses). Neglecting these feeble interactions and Z mixing, at low energies the HP interacts with the standard model particles only though kinetic mixing with the ordinary photon. The lagrangian density describing this low energy theory is where A µ , X µ are the photon and the HP fields, F µν , X µν their field strengths and j µ the ordinary electromagnetic current 2 , which for the purposes of this paper we can take as −eψ e γ µ ψ e (ψ e the electron field). The HP mass can arise either from the Stückelberg mechanism (see [38,42]) or from a hidden Higgs field developing a vacuum expectation value [71]. The latter case reduces to the former when the mass of the hidden Higgs field is much larger than the energies under consideration, but it is very different if the hidden Higgs is as light as the HP [50,71]. In this paper we focus on the HP phenomenology, which is equivalent to the Stückelberg case or the heavy hidden Higgs case. The A µ and X µ fields are non-orthogonal because of the kinetic mixing. The field redefinition diagonalises the kinetic part of the lagrangian, getting rid of the kinetic mixing, and reveals the interaction states A µ (photon-like, the eigenstate interacting with the electric charge) and S µ (sterile state that does not interact with the electric charge).

Photon "flavor" oscillations
The lagrangian (2.3) makes explicit that A, S are not propagation eigenstates in vacuum because of the non-diagonal mass term. The state radiated by electrons is pure photon-like (A) but, since it is not a propagation eigenstate, it will develop a S-component after some time and its magnitude is modulated by the time lapse. The picture is very much equivalent to neutrino oscillations: neutrinos produced in beta decay are produced in a flavour state (electron-flavor) and they oscillate into muon-flavor neutrinos, for example. In complete analogy with 2-flavor neutrino oscillations, the probability of finding a HP after a distance L from the photon production region is given by where ω is the photon/HP energy and we have assumed that HP are relativistic ω m. The role of the mixing angle in neutrino oscillations is played in this case by the kinetic mixing, χ. Since this has been already constrained to be very small, see Fig. 1, the probabilities of A → S transitions are minute and we always work at first order in χ. Since the difference between the states X and S is of order χ, we don't need to distinguish between them. When we compute oscillation probabilities, we refer to both simply as HPs.
The situation changes in a medium because photons (the A state) refract and get absorbed. The refraction and absorption properties in a linear medium can be casted in the form of an effective photon mass, which changes the definition of propagation states and thus alters the oscillation probability. This is again analogous to the case of neutrino oscillations, although in this context one speaks of "matter potential" instead of effective mass.
The probability of a photon of frequency ω to oscillate into a HP in a homogeneous absorbing medium has been computed in a number of papers in different frameworks: using the equations of motion in [49], Feyman diagrams and thermal propagators in [48] and thermal field theory in [46]. All these calculations lead to the same result in the small mixing regime, where the complex polarisation tensor Π = Π r + iΠ i (photon self-energy in the medium) depends on ω and can be understood as an effective photon mass m 2 γ and an absorption coefficient Γ which account for refraction and absorption (together with stimulated emission) of photons in the medium under consideration. Alternatively it can be expressed as a function of the complex index of refraction, N = n − iκ, Π = ω 2 (1 − N 2 ) 2ω 2 (1 − n) + i2ω 2 κ.

Inhomogeneous medium
In an inhomogeneous plasma, the probability amplitude of γ →HP conversion after a length L can be written as an integral over the putative photon trajectory r = r(l) as Here, ϕ is the the phase difference between the HP and photon waves integrated along the line of sight, in short, the number of γ →HP oscillations, and τ is essentially the optical depth. The physical interpretation of this formula is described in some detail in [19]. It is based on the perturbative expansion for photon-axion oscillations derived by Raffelt and Stodolsky [72] but applied onto the photon-HP system that was presented in [46] (kinetic approach of Sec. 3). At first order in χ, photons and HPs are propagation-eigenstates which can convert into each other with a a probability (amplitude) per unit length given by χm 2 /2ω. The transition can happen after any length l between the photon source and the (hypothetical) HP detector. The integral over the path-length l reflects this fact. The e iϕ(l)/2 factor is the phase difference between the photon and HP waves accumulated up to the length l.
Conversions after different lengths can interfere constructively or destructively. The factor e −τ (l)/2 reflects the absorption (and stimulated emission) of the photonic wave before reaching the conversion point at a distance l. See [19] for more details. From (2.7) and (2.8), it is straightforward to derive the two formulas for the probability in homogeneous media shown before. In vacuum, we have m 2 γ , Γ → 0 and thus while for a homogenous medium (2.10) This formula has to be a good approximation to media in which m γ , Γ change very little in an absorption length. The general formula (2.7) allows us to evaluate its limit of validity. Let us then consider a linearly changing medium, Expanding (2.7) up to linear order in the derivative we find Defining the mean-free-path as λ = Γ −1 0 and an oscillation length as λ osc = ω/(m 2 γ −m 2 ) the criterium reads approximately In this case ϕ(l) = , so the condition requires that ϕ is a linear function of l during a mean-free-path or an oscillation length, whatever it is smaller, up to corrections smaller than 1. The same of course, has to be satisfied by the optical depth function τ .
As it turns out, the condition (2.14) is satisfied almost in all the conditions of our interest because the typical length scale of an oscillation in the parameter range of our interest is extremely small, due to the large values of m 2 γ in the solar interior. Even the vacuum oscillation length, λ osc,vac ∼ 0.01 km ω 2 eV is much smaller than the characteristic density and temperature scale heights in the Sun, hundreds of kilometers at the surface and much longer in the solar interior.

Optically thick
The only obvious exception to the validity of (2.5) happens in regions where m 2 γ m 2 to good accuracy and the oscillation probability (2.5) is resonant. Close to a resonance, λ osc grows very large and eventually larger than the mean-free-path (which is typically longer than λ osc ). In the deep Sun, λ is extremely small and thus (2.14) is nevertheless satisfied. Resonances for which (2.14) holds can be called optically thick because their typical extent is and therefore (2.14) implies Γ∆r 1. Thus, (2.5) is an excellent approximation for the solar conditions except for what we shall call optically-thin resonant regions (Π 0 λ 2 /ω 1). Fortunately, for these cases we can also also derive a simple analytical formula.

Optically thin
The mean free path λ grows increasingly large as we exit the Sun, so the criterion (2.13) is eventually violated. The effect is exacerbated by the fact that near the solar surface, the density decreases much faster than in the interior and therefore so does Π .
Consider (2.7) and focus on a photon trajectory r(l) which intersects a resonance, i.e. a point where m 2 γ (r(l s )) ≡ m 2 γ (l s ) = m 2 . In this point dϕ/dl = ϕ = 0 and the integral has a saddle point, which typically gives the dominant contribution to the probability amplitude. A classical saddle point approximation gives where prime denotes differentiation with respect to l, ∆ ≡ l s |ϕ | and C(∆) is a smooth step-like function ∼ Θ(∆), which can be better approximated as The conversion probability can then be estimated as (2.20) The interpretation of this formula is very interesting. Most of the oscillations cancel out in the integral but there is a region of size (δl) 2 ∼ (ϕ ) −1 around l s where the phase has almost a constant value e iϕ(ls) and the amplitude has decreased by a factor e −τ (ls)/2 . The saddle point approximation singles out this region around the resonance as the most important contribution to the amplitude. In the particle conversion language, the photon→HP conversions along the path interfere destructively due to the fast varying (but smooth) function ϕ(l), except those happening in the region δl, which interfere constructively and thus dominate the probability.
Formula (2.20) is extremely simple and useful for our purposes. For the sake of completeness let us briefly discuss its limitations. The formula implicitly assumes that the saddle point dominates the integral. However, it is clear that if we separate enough from the resonance, at some point the exponential suppression due to the optical depth τ (l s ) will suppress the contribution from the saddle point below the contribution from the first mean-free-paths. In other words, for sufficiently large l s we are so far away from the resonance that λ osc is again small and the probability should again tend to formula (2.5) because (2.20) is exponentially suppressed at large distances. We expect a smooth transition between the two formulas. The situation is exemplified in Fig. 2 where we show examples of thin and thick resonance regions in the Sun. The left figure shows an example of photon→HP conversion probabilities around a thin resonance region. In this example, the photons originate at the given depth and move leftwards towards the solar surface. The red line shows a realistic approximation of the full result. In the deep Sun and far from the resonance it coincides with the homogenous approximation (2.5) (shown in blue) but as the production point gets closer, the contribution from the resonance region starts to dominate and grows exponentially (here the mean-free-path around the resonance is λ ∼ 300 m). Our probability formula (2.20) captures this exponential growth. Once the resonance is past, the probability returns quickly to the homogeneous result. The oscillations seen at the beginning of the exponential growth originate from interference of the local and saddle point contributions and they are never relevant as they tend to average out when frequency averaged or when the probability is integrated over some region of the Sun as when we compute the HP flux of the entire Sun.
In contrast, the right figure shows a thick resonance region where the criterion (2.14) is satisfied and photons oscillate into HPs before realising any inhomogeneity in their index of refraction or absorption (at the resonance centre λ ∼ 3 m).
Note finally that in order to use (2.20) for ∆ < 0 (l s < 0), i.e. for photons moving away from the resonance, we have to taylor a small modification. The required change is to use τ (l s ) = 0 for l s < 0 because the region expected to contribute to the integral is the closest to the photon emission, thus τ (l s ) → τ (l s )Θ(l s ). As ∆ becomes negative the photon path contains less and less of the maximal mixing region and the amplitude decreases as ∝ 1/∆. In this regime we can recover explicitly (2.5) in yet a different fashion. Aside from irrelevant phases we find where for the last expression we have used the Taylor expansion of m 2 γ around the resonance m 2 γ (l) = m 2 + m 2 γ (l s )(l − l s ) + .... We obtain this expression because in this limit, most of the conversion probability comes from the first oscillation, as in the homogeneous case.

General aspects
The emission rate of transversely polarised HPs of energy ω per unit volume at a given position, r 0 , of an inhomogenous plasma in local thermal equilibrium (LTE) can be written as the photon production rate Γ P times the conversion probability P (γ → HP) [49] integrated over the different directions in which a photon can be produced, Note that the conversion probability P = P (γ → HP) depends in principle upon the whole trajectory r(l) ∼ r 0 + lk and in particular on the direction of the momentum k. The photon production rate only depends on the creation point, as follows from our assumption of LTE. The factor of 2 accounts for the two transverse polarisations. In LTE, the photon production rate is related to the absorption rate Γ A by detailed balance, Γ P = e −ω/T Γ A , and the imaginary part of the self energy, Γ = Γ A − Γ P , so that with T = T (r) being a smoothly-varying plasma temperature. The specific HP flux on Earth is the integral of the volume emission of the Sun divided by the surface of a sphere of radius the Sun-Earth distance, R Earth 149.6 × 10 6 km, where we have used the HP dispersion relation ω 2 = |k| 2 + m 2 .
Let us now assume the Sun to be spherically symmetric and postpone the discussion on inhomogeneities. Then the probability P only depends on the radial position, r, and the azimuthal angle, θ, which defines the photon/HP trajectory, see Fig. 3. Integrating over the remaining angular variables we find In order to compute this integral, we need to know the functions Γ(ω, r), T (r) and P (ω, r, θ), for which we also need m 2 γ (ω, r). In Sec. 3 we describe how to build these expressions from a solar model and in Sec. 4 we perform the relevant computations to obtain the solar HP flux and discuss the results.

Resonance domination
It turns out that for values of ω in the visible range, the conversion probability is very peaked around a spherical shell located at radius r * defined by where the photon-HP oscillations are resonantly enhanced. The emission from this region can easily dominate the emission from the rest of the Sun, making possible a humongous simplification of the calculations. Indeed we can perform analytically the remaining two integrals of (2.26) to obtain a surprisingly simple result. The derivation will take the rest of this section.
We have simple probability formulas which apply when the resonance region is optically thick or optically thin. The criterion we shall use to identify them is Note that the gradient depends on the azimuth of the trajectory and thus resonances are thin or thick depending on cos θ. It is interesting thus to focus first on the r-integral, which we do in the following.

Optically-thick resonance
If the resonance is optically thick, the γ →HP conversion probability is well approximated by (2.5) and the r-integral is which actually does not depend on cos θ itself. In the resonance, the probability is enhanced by a factor m 4 /(ωΓ) 2 with respect to the vacuum case and much more if we compare it with very dense regions where m 2 γ (r) m 2 . Since usually we have m 2 γ ωΓ, the integral so strongly peaked that we can approximate the Lorenzian shape of the probability by a Dirac delta to obtain (2.32) The explicit dependency on Γ drops out [46,48,49] and the resulting flux only depends upon our solar modelling through m 2 γ (ω, r) and T (r * ).

Optically-thin resonance
If the resonance is optically thin we can use (2.20) for the r-integral, which now depends explicitly on θ because of the cos θ factor in (2.30) and the optical depth to the resonance, The r-integral is where we have taken C(∆) = Θ(∆), which we write as Θ ≡ Θ ((r * − r) cos θ), ensuring that only photon trajectories that cross the resonance are included. If the photon is produced at r < r * only trajectories with cos θ > 0 will cross the resonance, the opposite case being when r > r * . The integral is dominated by the region where the optical depth to the resonance is O(1), typically a few mean-fee-paths. Assuming that r and T do not change much in that narrow region we can substitute for the values at the resonance, r * , T * and perform the integral explicitly where τ * is the optical depth of the resonance with respect to the surface. If the resonance lies deep enough inside the Sun (τ * ∼ 3 − 4 or so), HPs converted from photons travelling through it from the inside or from the outside contribute an equal amount to the total HP flux. If it lies close or in the photosphere, the amount of photons produced outside and traveling inwards decreases very much and we are left only with the HP converted from outgoing photons. Remarkably, again the resulting flux does not depend explicitly on Γ or other properties of the solar model other than T and m 2 γ . We have then Remarkably, it is given by exactly the same expression than for the optically thick resonance (i.e. as long as τ * 3 − 4).

Master formula for resonance emission
If we assume that the resonance region is either thick or thin or, more specifically, that the values of cos θ for which the resonance is neither are not quantitative relevant for the cos θ integral, the cos θ integral is trivial where δθ = thin d cos θΘ(− cos θ) corrects for the HPs produced from photons originated at r > r * which travel inwards the Sun. The 1 in the formula would be symmetric result, which is to be corrected when e −τ * is not negligible. In practice, δθ 1 for all practical purposes. Our final formula for the HP emission from the entire resonance region is The formula was already found in [46,48,49] under the assumption that the resonance region is optically thick (extremely good approximation in the solar interior, where these references were focused). In this paper we have shown explicitly that its validity extends to optically thin resonances and for the mixed case where the resonance region is optically thin for cos θ ∼ ±1 and thick for cos θ ∼ 0. For that, we needed to assume that the regions in θ for which the resonances are neither thick or thin do not change the behaviour. I have performed many explicit calculations and cross-checks that show that this is the case, but at the moment I have no general proof. However, in the simple case of constant r, T, Γ, dm 2 γ /dr in the resonance region (which is not a bad approximation in the Sun) one can perform analytically all the integrals and prove that there is no funny behaviour when the resonance is neither thin or thick. The calculation is shown in appendix A.

Angular distribution of the signal
The specific HP flux at earth (HPs per unit area, time, energy and stereoradian) along a line of sight intercepting the Sun is where r = r(s) with s the line-of-sight distance from the Earth to the production point (dashed line in Fig. 3). Defining r min = R Earth sin ψ as the impact parameter of the trajectory we change the integration variable to the radial coordinate r If the trajectory intercepts the resonance region, the bulk of the emission can be readily evaluated by similar calculations which took us to (2.40). The key simplification is to use the optically thick resonance formula (2.5) for the probability because either r 2 − r 2 min is relatively constant during the integral (and then optically thin and thick resonances give the same result) or r * ∼ r min and then cos θ ∼ 0 and (2.5) is justified again. Using r ∼ r * in the numerator and r 2 − r min ∼ (r * + r min )(r − r min ) in the denominator, the integral can be readily evaluated as dΦ(ω, ψ) dωdΩ where the angular distribution is given by where the angle at which ψ is tangential to the resonance shell is ψ * r * /R Earth and the width, which we have assumed to be small, is given by Moving from the solar center outwards, the angular distribution grows as ψ → ψ 0 and peaks at ψ * − ψ ∆ψ * , where the line-of-sight is tangential to the resonance shell and thus benefits from more resonant emission, and then decrease extremely fast. The peak has a 1/ √ ψ − ψ * behaviour and of course does not dominate the integral.
Since the resonance can be extremely sharp it might turn out that cannot be resolved by a telescope. In this case, the resonance will be broaden by the finite resolution of the apparatus. The angular integral of the resonance is independent of the resonance width as long as it is very small (2.47)

Basics
In order to compute the solar HP emission, we need to model the refraction and absorption properties of light in the solar plasma. For the latter, very exhaustive studies exist since absorption determines the radiative energy transfer inside of the Sun, which in turn determines the solar structure. We highlight the Opacity Project [73-76] and Los Alamos opacity code LEDCOP opacities [77][78][79]. To the best of our knowledge, no study of the refractive index throughout the solar plasma exists in the literature. In principle, the real part of the index of refraction can be obtained from the imaginary part with the Kramers-Kronig relation. Unfortunately, the existing data is often smoothed over frequencies and it is only available for a small number of density and temperature points. A rough interpolation would introduce large errors in the determination of the resonance region so we have to follow a different procedure.
We can calculate explicitly the most relevant contributions to refraction and absorption as a function of temperature, density and composition so that we have a very smooth map of them inside of the Sun (as smooth as the solar model we might use). These are the contributions from electrons either free or bound in H atoms. The effects of electrons bound in Helium and heavier atoms are subdominant for refraction, but can be relevant for absorption around frequencies corresponding to the strongest atomic transitions. This contribution is then taken from existing opacity calculations to avoid the extremely involved atomic physics. Indeed, we can use OP opacities tables, conveniently provided for each metal separately and thus allowing arbitrary mixtures.
The solar plasma consists on electrons (free or bound in ions) and nuclei, which being much heavier interact weaker with light than the former and can be neglected. Unbound (free) electrons contribute to the polarisation tensor with a simple term n free e (3.1) whose real and imaginary parts we recognise as the plasma frequency squared and the absorption coefficient due to Thomson scattering. Electrons bound in H atoms can be in first approximation modelled by a set of oscillators whose contribution of the polarisation tensor is where n H 0 is the number density of neutral H atoms, Z n the probability of finding the bound electron with principal quantum number n and the last sum is over resonant transitions n → n , which happen at frequencies ω r = Ry 1 The oscillator strength f nn is summed over final states and orbital quantum numbers l and l = l ± 1 for electric dipole transitions. We neglect fine-structure corrections as we are not interested in fine details of the spectrum. The width of the resonance is given by the natural line width 2αωω r /3m e . Impact broadening due to proton collisions through the Stark effect is taken into account in the quasi-static approximation by folding each resonance contribution with a Holtsmark distribution following [80]. We don't consider collisional broadening from electron collisions nor Doppler broadening, which are not relevant at the level of accuracy required. Broadening is relevant close to the resonant transitions and in the wings of absorption lines for Γ, but not for m 2 γ . Note that free electrons behave as an additional bounded species with f = 1 and ω r = 0.
The above expressions account for the most important contributions to refraction but are not enough for absorption. The leading contribution to photon opacity in the deep Sun comes from photon absorption during electron-proton scattering, γ + e − + p + → e − + p + . In outer shells, also the photoelectric effect, γ + H * → e − + p + , is important. Usually these reactions are known as free-free and bound-free processes, respectively. Their contribution to Π i can be casted in the following compact form where F ff , F bf are thermally averaged Gaunt-like factors, introduced as corrections to the classical result which depend mildly on frequency and atomic details and can be found for instance in [81]. Close to the threshold F bf ∼ 1 so we neglect it. As for F ff , for our purposes it is sufficient to consider the Born-Elwert approximation [82] with a simple screening prescription 2m e T and we take k D to be the Debye-screening scale given by k 2 D = 4πα α Q 2 α n α /T where the sum extends to all charged particles (mostly electrons, protons and some He atoms). In the bound-free expression, E n = Ry/n 2 is the energy of the n-th energy level.
Despite its low density, the negative H ion, H − , plays an important role in the opacity near the photosphere at near IR and visible frequencies [83]. The photoionisation cross section of H − does not have a sharp edge at its ionisation threshold ω = E − = 0.75 eV because the electrons are ejected in a p-wave [83,84] due to the neutrality of the H 0 final state. Instead of the ∝ Θ(ω − E n )/ω 3 behaviour of other contributions to Π bf one gets ∝ (ω − E − ) 3/2 /ω 3 . We have taken the photonization cross section from [85].
These contributions have associated refractive parts, which can be computed through the Kramers-Kronig relations, neglecting the ω dependence of the Gaunt factors. The boundfree contribution of H 0 turns out to be of similar size to the bound-bound contributions and has to be included. Due to the low density and the smoothness of the threshold the refractive part associated with the photionization of H − turns out to be negligible. Neglecting the frequency dependence of F bf we find Finally, we want to include the absorption coefficient due to metals. In principle, suitable generalisation of the above formulas is possible and regularly performed for opacity calculations. For the scope of this paper it is enough to use monochromatic opacities present in the literature. We chose the detailed opacities of the Opacity Project from [76], which are available in tables for different temperatures and electron densities that can be interpolated at will. The absorption coefficient is tabulated for each element as an effective cross section that includes the contribution from scattering, bound-bound, bound-free and free-free reactions involving the element Z, σ Z (ω) = Γ Z /n Z (1 − e −ω/T ) where n Z is the density of the corresponding nucleus. We thus build (3.8) The polarisation tensor Π is the sum of all the above contributions,

A model for refraction and absorption in the Sun
The mass density ρ, chemical composition and temperature T as a function of solar radius can be taken from a standard solar model calculation. The state of the art on solar modelling is able to fit all available solar data from helioseismology and neutrino flux detectors with typical percent accuracy. Small discrepancies with observations still exist [86] after the recent revision of solar abundances of CNO species [87], which lowered slightly the opacity, but they do not compromise the degree of precision of our calculation. Interestingly, the authors of [88] proved that the existence of hidden photons cannot possibly influence the determination of the light elements in the Sun.
In this paper we make use of the Saclay seismic solar model [89,90] available at [91] which has excellent detail in the outer layers of the Sun, our primary concern. The mass density and temperature profiles are shown in Fig. 4. For future reference, Fig. 5 shows the outer layers of the Sun in more detail as a function of the depth towards the solar core.
The surface chemical composition from [87] is shown in Fig. 6. Hydrogen is ∼ 10 times more abundant than helium, ∼ 1000 times more abundant than CNO elements and more than 10 4 than higher−Z elements. The composition is relatively constant in the solar interior but we take into account diffusion from the solar model AGSS09 [92], available at [93] (the solar model [89] does not provide the chemical abundance as a function of radius).
The most relevant atomic species is of course, hydrogen. As we will see further on, the resonance region for low mass HPs happens close to the photosphere but still in the solar interior. In this region the temperature is so low that He atoms are not ionised (their first ionisation energy is ∼ 25 eV) and their contribution to free-free transitions is negligible. Also, the main atomic transitions and ionisation thresholds lie in the far UV, and therefore do not affect visible light absorption. For the same reason, their contribution to refraction is also   negligible. Note from (3.2) that the contribution to Π r of far UV resonant frequencies in the visible is ∝ −n Z (ω/ω r ) 2 and therefore those of helium are much less important than those of hydrogen. Higher Z elements (often referred to as metals) are too scarce to compete with hydrogen in absorption in free-free collisions but contribute to bound-bound and bound-free absorption because they have resonant transitions in the visible and ionisation thresholds in the not-so-far-UV. The effects, if competitive with hydrogen have to be necessarily localised in frequency to very narrow intervals around the resonances and thresholds. Since we are interested in the continuum flux rather than a spectroscopically resolved spectrum we will neglect these contributions. The only effect in which metals have a leading role is to set the free electron density in the photospheric layers of the Sun where hydrogen is highly neutral. Only in this restricted sense we include metals in the solar refractive model.

Partition functions and atomic ionisation
The density of free electrons and of the 1s state of neutral hydrogen are the most influential factors that set the value of m 2 γ . Through this dependence, they determine the locus of the photon-HP resonant conversion region (and hence the temperature T of the emission), and also the size of the conversion region, which sets the absolute value of the flux. Thus, for our purposes it is of capital importance to have a robust estimate of these quantities. The occupation probabilities, Z n , and the neutral H fraction cannot be accurately estimated with the partition function of an ideal gas and the Saha ionization equation. The first requires an ad hoc cutt-off to the number of excited states, and we do not know a-priori if this is justified. Excited states could be very populated at large temperatures and contribute to m 2 γ as much as the 1s. The second is well known to fail at large densities and not too large temperatures, predicting too little ionisation, even in the solar centre!. Fortunately, the calculation of the ionisation equilibrium is a fundamental ingredient of the equation of state of astrophysical plasmas and radiative opacity, so a large body of literature exists on the matter. The problems of the divergence of the partition function and the Saha equation are both cured by include interactions between the plasma species, which alter the bound states when the density gets closer to atomic density. Here, we use the partition function of Hummer and Mihalas (HM) [94], where the factors w n are called "occupation probabilities" and encode the non-ideality of the gas. The most relevant effect to be included are the perturbations on the atomic states by the slow-varying electric fields of the plasma ions (protons in our case). Bound states suffer Stark shifts in the binding energies which drive ionisation very efficiently. In particular, there is a certain critical field above which a given bound state is has a lifetime smaller than its orbit period and it effectively destroyed. These effective occupation probabilities thus measure the probability of the atoms to be in a perturbing field such that a given bound state does still exist. Using a Holtsmark probability distribution for the magnitude of the fields, they find where Q(x) is the cumulative of the Holtsmark distribution, n p is the proton density and K n = 16n 2 (n + 7/6)/3(n + 1) 2 (n 2 + n + 1/2) for n ≥ 3 and K n = 1 for n ≤ 3. The Saclay solar model [89] uses the OPAL equation of state (EOS) of Rogers and Iglesias [95], which is not based on the HM chemical picture but on the physical picture (see [96] for the derivation of occupation probabilities). The Opacity Project EOS is instead based on the much simpler HM picture and it has been shown to agree well with the OPAL [97]. Therefore, by using the HM formalism we simplify much our calculations without the danger of large inconsistencies on the chosen solar model. Note finally, that the chemical picture of HM was designed in principle to be valid at the small densities relevant in stellar envelopes, ρ 10 −2 g/cm 3 but it was found to be valid up to much denser regions, and in particular, the whole Sun [98]. An update of the micro-field distribution including particle correlations showed an excellent agreement with OPAL up to the solar centre [99], but for our purposes we do not need such levels of precision.
Following HM, the ionisation equilibrium is derived from the minimisation of the Helmholtz free-energy where every atomic state is treated as a different species. We consider the states of the hydrogen and helium and the first ionisation of metals. There are three states of H (H − , H 0 and H + ≡ p, i.e protons), three of He (He 0 , He + and He ++ ) and two for each metal. The effects of H bound states are only important in the surface, where excited states are scarce and therefore their correlations with the proton density irrelevant for the free energy. In this lucky situation, the minimisation of the free-energy, F , subject to the stoichiometric relations where we have defined the usual partition function normalised to the ground state energȳ Z = e E 1s /TZ (recall that E n > 0 in our convention) and the Debye correction involves the Debye screening scale defined before. The equations for H − , He and metal ionisation are completely analogous. The 2 Saha equations of H, the 2 of He and the one for each metal are to be solved self-consistently with the constraints given by the solar abundances (n H,total = n H + + n H 0 + n H − for H and n He,total = n He ++ + n He + + n He 0 for He) and the expression of the free electron density n free e = n H + − n H − + 2n He ++ + n He + + Z n Z + . (3.13) The resolution of this system outputs the required ion densities and bound-state probabilities Z n , which we show in Fig. 7.
Outside the photosphere, neutral H and He are the most abundant species. The free electron density has dropped much, although not as much as the proton density due to The ionisation of H is around 90% around ρ ∼ 10 −3 g/cm −3 and that of helium around ρ ∼ 10 −2 g/cm −3 . Let us warn that the free electron density profile obtained by this whole procedure does not fit well with the one provided by the solar model [89] below ρ ∼ 10 −6 g/cm 3 or so. Our calculations show a decrease of two orders of magnitude around ρ = 3 × 10 −7 g/cm 3 , which is much milder in the solar model table provided in [91]. The reason is most likely that the available model used in [89] employs a basic Hopf atmosphere model, see [90]. The results of our calculation agree perfectly with the solar atmosphere model of Kurucz [103], which the authors of the Saclay solar models employ in the later publication [90] to improve the agreement with Helioseismological data. This constitutes a reassuring test of our calculations.

Results
We can now construct the functions m 2 γ (ω, r) and Γ(ω, r) that we require for computing the solar HP flux. We show some relevant plots of the dependence of m 2 γ (ω, t) on r and ω in Fig.  8. In the upper plots we see the frequency dependence of the different contributions. Free electrons provide a frequency independent contribution (blue line) while atomic resonances at ω r give positive contributions for ω > ω r and negative for ω < ω r . Negative values are shown as dashed lines in the log-plot. The most notable line is the Ly-α at ω ∼ 10.2 eV although all the Lyman series has visible contributions. The value of m 2 γ is determined thus not only by the density but very importantly by the ionisation fraction. In the surface, i.e. just above the fast drop of density visible in Fig. 7 (ρ ∼ 3 × 10 −7 g/cm 3 ) most of the H is neutral and m 2 γ is negative in all the visible spectrum. Moving deeper in, at a density ρ ∼ 5 × 10 −7 g/cm 3 the ionisation fraction is already 10%. Since the neutral H contribution is ∝ −ω 2 it becomes ineffective at low energies and such a small free electron density is able to make m 2 γ positive in the red part of the spectrum. As we move inside the Sun, the ionisation fraction increases and neutral H becomes scarce. The region of negative m 2 γ gets displaced deeper into the UV. At ρ = 8 × 10 −5 g/cm 3 it is only the ∼ 9 − 10.2 eV range. These trends can be seen in the lower plots of Fig. 8 where the dependence with the solar density is shown for a couple of frequencies. The positive free electron contribution drops suddenly at the photosphere with the ionisation fraction while the negative contribution from neutral H decreases softer and eventually dominates. Note that for these frequencies there is always a point in the Sun where m 2 γ = 0. Since the negative contribution is frequency dependent this point depends on frequency. This point will be very important for the discussion on photon-HP oscillation resonances. Note that as we consider smaller ω the neutral H contribution (red-dashed line) decreases as ω 2 and the point when it crosses the free electron (blue), i.e. the point where m 2 γ = 0 moves out of the Sun. Higher energies have m 2 γ = 0 deeper in the solar interior because of their larger neutral H contribution. Interestingly, due to the sharp drop of the free electron density at ρ ∼ 3 × 10 −7 g/cm 3 a sizeable range of frequencies will have m 2 γ = 0 around that region. The imaginary part of the self energy (the absorption coefficient up to a small correction) is shown in a few plots in Fig. 9. In general it agrees well in the range of interest with the calculations of the OP interpolations and crosschecks with the LEDCOP database [77][78][79] at the low densities relevant for this work.

Solar hidden photon flux: Numerical results
The model for refraction in the solar interior build in the last section allows us to compute the solar flux of HPs by direct integration of formula (2.7). This is a very time-consuming operation, which fortunately it is not necessary because the simple formula (2.5) is valid for most of the Sun, with the only exception of optically thin resonance regions where (2.20) can be used. The solar HP flux in the visible is dominated by resonant production, which we discuss in detail next. Later, we integrate (2.5) in the whole mass and energy range to compute the non-resonant contribution and present our atlas of solar HP emission. Finally, we estimate the corrections to the 1D solar models used by computing the emission from a 3D time-dependent solar atmosphere model.

Resonance region contribution
In section 2.5.2 we showed that the formula (2.40) provides an accurate description of HP flux from the resonance region. In order to evaluate this expression we need the position of the resonance as a function of the frequency and the HP mass, r * = r * (ω, m). We have solved numerically the equation m 2 = m 2 γ (ω, r * ) and present in Fig. 10 our results. In the regions shown, the resonance moves to low densities with decreasing HP mass and photon/HP energy. The first trend is obvious from the equation m 2 = m 2 γ (ω, r * ) because m 2 γ is proportional to the densities of charged particles. The second follows from the fact that, at low densities, m 2 γ has a positive contribution from free electrons and a negative one from neutral H which decreases with ω 2 . A decrease in the negative contribution has to be balanced by a decrease in the positive one, i.e. displacement towards lower densities where the free electrons become much more scarce.
Our results are shown in Fig. 11. In the frequency range shown, the flux peaks at low frequencies and decreases strongly with the mass (note that the fluxes in Fig. 11 are divided by m 4 for the sake of illustration). Let us first comment on the spectral shape. It is the result of the convolution of the factor √ ω 2 − m 2 /(e ω/T (r * ) − 1) and the derivative |dm 2 γ /dr| −1 which gives the typical size of the resonant region. Both are larger at low energies. The first term is relatively flat because, at the masses and frequencies shown, T (r * ) ∼ O(1) eV. Taking ω m we have √ ω 2 − m 2 /(e ω/T (r * ) − 1) ∼ T (r * ). The exponentially suppressed regime only starts to be felt at the highest energies and the threshold at the largest masses. The second term can be understood if we write a very schematic model for m 2 γ at low energies as where n t is the total number density of H atoms, X e = n free e /n t is the ionisation fraction and ω 0 ∼ 10.2 eV is the resonant frequency of the Ly−α transition. The first term is due to free electrons and the second to neutral H. The derivative with respect to the radius can be written as n t is a relatively smooth decreasing function of the solar radius and X e is essentially = 1 flat in the interior and drops exponentially fast near the surface (see Fig. 7 right). Resonances happening in the deep Sun have d log X e /d log r ∼ 0 and thus dm 2 γ /dr independent of ω. But for those happening near the surface, the d log X e /d log r term dominates. In that region, the resonance condition m 2 γ = m 2 is realised despite a large value of 4παn t /m e > m 2 by some degree of cancellation in the term X e − (1 − X e )ω 2 /ω 2 0 so that X e ∼ ω 2 /ω 2 0 and the derivative becomes suppressed at low energies. The suppression is not proportional to ω 2 because d log X e /d log r is very sensitive to r * , which increases for lower ω but, nevertheless, the general trend of smaller dm 2 γ /dr and thus bigger fluxes remains. Let us know shed some light on the dependence with the mass. At high masses, the resonances happen in the deep Sun and, according to (4.3), the derivative dm 2 γ /dr is proportional to m 2 γ m 2 so the flux should go as m 2 |d log n t /d log r| −1 , which when properly evaluated at the resonance point gives something ∝ m 3 . This trend shows in the mass region 0.01 − 0.1 eV in Fig. 11 (right) and it is only stopped at higher masses because we have focused in the ω ∼ O(eV) region and the kinematic threshold cuts the curves. At masses below 0.01 eV the resonances approach increasingly the photosphere and they cannot go further. This is because, for all the frequencies shown there is a point where m 2 γ becomes 0. HP masses of arbitrarily small mass, will have their resonance arbitrarily close to that region, but not further. In that limit, dm 2 γ /dr will be the same for all small masses (although it depends on the energy, as we have already explained). This explains the flattening of the low mass flux below m ∼ 0.01 eV in Fig. 11 (right).

The infrared rise
If we consider sufficiently low energies and low masses the resonant region starts to move out of the photosphere out into the open atmosphere. In Fig. 10 this is seen to happen in the region m 10 −4 eV and ω < 0.4 eV. In this region, the photon mass m 2 γ decreases smoother than in the surface, (see the electron density at ρ ∼ 10 −8 − 10 −7 g/cm 3 displayed in Fig. 7 (right)) and consequently the resonance region is larger (|dm 2 γ /dr| smaller) and the photon→HP probability gets boosted. For m = 10 −4 eV, the flux at ω ∼ 0.2 eV seems to be more than two orders of magnitude stronger than at 2 eV. But it is not clear that we can claim these results to be true because our solar model for refraction is probably not very accurate here. There are a number of effects that can alter our calculation. Let us discuss some of them. First of all, where the rise starts depends very much on the precise determination of the profile of the free electron density. In the atmosphere, the contribution of metallic donors is crucial, and we have performed only a rough estimate. Second, even a small number of non-H atoms with resonances in the visible or IR could in principle contribute more than the UV resonances of H to the index of refraction here. These contributions are negative and push the resonances again deep into the solar interior where the gradient |dm 2 γ /dr| is larger and thus the flux smaller. A simple estimate tells us that this uncertainty grows very much at low energies. Imagine that every metal, but not He, contributes one electron with a resonant transition at frequency ω m , then in the far infrared where m 2 γ,H is the estimation of the H contribution through the Lyman series. For ω m in the visible, this effect is negligible but if there is considerable structure below 0.3 eV the effect would be similar to the H. This rough estimate, shown with the aim of highlighting the typical orders of magnitude is enough to rise a voice of warning towards the refraction model validity in the IR. In the atmosphere, the density of H molecules like H 2 , CH, OH, is extremely small but they have absorption lines in the IR that one should in principle take into account. Another assumption which is prone to fail in the atmosphere is LTE, which we have used to compute all the abundancies and the radiation temperature. Finally, we will see that the spherically symmetric model fails to some extent to reproduce the true 3D behaviour of the solar surface. Due to these issues, in this paper we can make no claim of any rigour in the derivation of the IR flux produced in the solar atmosphere. Although it is clear that the typical flux grows if resonances happen in the shallower density profile of the atmosphere, the low energy limit of Fig. 11 have to be understood as order of magnitude estimates.

Spectral lines and the UV region
The results shown in Fig. 11 show only tiny spectral features, which prove that the role of atomic resonances in the visible is not crucial. Near a strong spectral line, the pattern of the flux is expected to change because, just below, the effective mass m 2 γ is more negative and above, more positive. Resonances are displaced towards the solar interior in the red part and towards the surface in the blue and consequently, the red part is more luminous that the blue part. An example, corresponding to the H-α line is shown 3 in Fig. 12. Even around this important line, the effects are at the 10 percent level. We will see later that O(1) fluctuations in the position of resonances close to the photosphere and their associated fluxes arise due to convection so these effects are completely unobservable. One can also see that the spectral feature softens as the HP mass grows and the resonance region enters deeper into the Sun, where the line is broader.
The behaviour around other lines is qualitatively similar to H-α. Lines of the Lyman series give the most spectacular effects in the UV, shown in Fig. 13. The Ly-α line at 10.2 eV is so strong that pushes the low mass HP resonances up to 5 Mm deep inside the Sun   at ω ∼ 10 eV and expels them outside of the Sun in the 10.2 − 11 eV range. The other Lyman lines have similar effects, bring the resonances some Mm inside of the Sun at their red sides and expel them to the solar surface in the blue sides. The resonant flux in the UV varies very violently with energy, dropping by many orders of magnitude in the blue sides of lines. Indeed, in these dips the resonance flux becomes subdominant with respect to the non-resonance flux as we shall see in section 4.2. The situation softens progressively for increasing HP masses and it is very smooth above m ∼ 0.4 eV.

Width of the resonance region
If the resonance is optically thick, it follows from (2.31) that the FWHM of the resonance is ∆r thick ∼ ωΓ(r * ) dm 2 γ dr −1 r * [49]. When it is optically thin, it is given by the mean-free-path around the saddle point, ∆l ∼ 1/Γ(r * ) but since l ∼ ∆r/ cos θ, ∆r thin ∼ cos θ/Γ(r * ) it is azimuth dependent. We can then estimate the resonance width as which is shown in Fig. 14 as a function of the HP energy. The red range around the resonance position represents ∆r thin for cos θ = ±1 and the turquoise region ∆r thick . In resonance regions close to the solar surface, scale heights are very short and thus ∆r thick decreases while mean-free-paths become longer. These resonance regions are optically thin, except near strong absorption lines. Note that in these resonances, trajectories with cos θ 0 are thick because traveling in a direction perpendicular to the radial direction, the photon/HP sees the same properties of the Sun before being absorbed. So, in reality, the thickness of resonances in the red regions vary from the red to the turquoise range depending on cos θ. The size of the saddle point region, δl ∼ |ϕ | −1/2 = dm 2 γ /dr2ω −1/2 , is smaller than 100 meters in the region shown. Deeper than about 0.4 Mm, all the resonances are optically thick. Comparison of Fig.  14 with the temperature and density profile of Fig. 5 shows that that in general, resonances are extremely thin with respect to the temperature scale height, even around absorption lines or phoionisation thresholds. Assuming constant temperature and density within the resonances seems a reasonable approximation in this energy range.

Non-resonant contribution (bulk of the Sun)
So far we have focused on the resonance region and its contribution to the HP flux. We have advanced that it is the most relevant contribution in the visible range of energies. In this section we justify why this is so and then compute the non-resonant flux in the UV and above where it can dominate.
Consider the integral over the Sun that gives the total HP flux on Earth for a given HP energy, (2.26). The locally homogeneous plasma approximation, (2.5), is an extremely good approximation for the photon/HP conversion probability, except for optically thin regions, which however are so thin (like any optically thick) that do not make a difference in this discussion. Moreover, we know that once integrated over the Sun, the result is independent of the optical thickness of the resonance. After performing the now trivial cos θ integral, the HP flux in this approximation is (4.6) There are three qualitatively different regions in this integral. The most notable is the resonance region, of typical size ∆r ∼ ωΓ|dm 2 γ /dr| −1 R Sun where the conversion probability is maximal ∼ χ 2 (m 2 /ωΓ) 2 because typically m 2 γ = m 2 ωΓ (with exceptions being strong absorption lines, for which the resonance is proportionally wider).
Outside the resonance region, m 2 γ m 2 , the conversion probability is like in vacuum ∼ χ 2 , the region is bigger than the resonance region but temperatures and production rates are very small so that we shall not expect competitive rates from here.
Inside the resonance region r r * we have m 2 γ m 2 and thus (m 2 γ ) 2 +(ωΓ) 2 m 4 and the conversion probability is suppressed with respect to the vacuum case ∼ O(χ 2 ). However, this region has the typical size of the whole volume of the Sun, and the factor Γ e ω/T −1 is larger because so are the temperature and the electron/proton densities that determine Γ. One could be tempted to think that, even if the probability is suppressed, the net flux is competitive with the resonance. What is stronger, the Γr 2 /(e ω/T − 1) enhancement or the (m 2 γ ) 2 + (ωΓ) 2 suppression? The situation at visible energies is depicted in Fig. 15 flux originating from the bulk Sun inside the resonance region drops dramatically below the resonance and then stabilises to a much smoother decreasing contribution. In the solar interior the refraction is dominated by free electrons (m 2 γ = 4παn free e /m e ) and the absorption by free-free transitions (Γ 64π 2 α 3 √ m e (1 − e −ω/T )n free e n p /(3m 2 e ω 3 √ 2πT ) neglecting F f f and screening) from which we can already see that the density dependences cancel out in Γ/m 4 γ (n free e ∼ n p ) and the resulting formula depends smoothly on the solar parameters. We where T, ω in this formula are in eV units. These estimates are shown as dashed lines in Fig.  15 and fit very well the full formula at depths larger than 10 Mm. With R Sun = 695 Mm, one could think that this contribution is comparable to the resonance contribution for the highest HP masses (O(eV) in this context of visible energies), however, the factor of √ T ω in the denominator suppresses enough the contribution.
Finally, we have performed numerical integrations of (4.6) in the 1−10 4 eV energy range in a generous range of masses to compare with the resonance estimate presented before. The atlas of HP solar emission is presented in Fig. 16. For each HP energy and mass we show the full emission from the integral (4.6) (solid line) and the contribution from the resonance region estimated by (2.40) (dashed line). The solid lines are coloured-coded according to the contribution that dominates the overall flux: black where the resonance accounts more than 60%, red when it is less than 40% and an interpolation in blue in between. The red colour reminds that these results depend on our model of Γ.
In the visible, the results agree notably well with the resonance contribution. Only above m > 0.1 eV the all Sun integration rises the prediction somehow, up to 40% in some cases. However, most of this extra flux probably comes from the tails of the resonance region, particularly the inner one where the temperature T and production rate ∝ Γ are larger and can rise somehow the flux. So this deviation has to be understood not as a contribution from the bulk of the Sun but as a correction to the thin-resonance condition that we have used.
In conclusion, the resonance region dominates the production of solar HPs in the visible and provides a precise prediction independent on the uncertainties on Γ. In the UV, the situation changes very much. For low mass HPs (below m ∼ 10 eV), if the resonance takes place, it does so in regions where the temperature is O(eV). The resonant production of this high energy HPs starts to be exponentially suppressed and thus the bulk contribution competes. In the near UV, below ω ∼10 eV the overwhelming number of atomic transitions of metals rises the resonance prediction. In the Lyman region (10-14 eV) we see the structure advanced in section 4.1.2. In the blue part of the Lyman atomic transitions, the resonance is not possible and the flux is lower and dominated by the bulk (red dips), while in the red part there is still some resonance emission, which turns out to be of the same order than the bulk emission (probably also because the resonances are not thin). In the far UV the exponential suppression is huge and we only see the bulk emission. For instance, we can see the 1s-2p transition line of singly ionised helium at ω ∼ 40 eV and guess a multitude of smaller lines from metals. From ω ∼ 300 eV we see strong line emission in the soft X-ray region that decreases at 3 keV but shows up to 6-8 keV with the strongest Iron lines.
As we consider larger HP mass, the transition between resonance domination at low energies and bulk domination at the highest displaces towards higher energies because the resonance region is located deeper into the Sun and happens at higher temperatures. For instance, HPs with masses around and above 10 eV have resonance regions where the temperature is high enough to dominate the full emission up to ω ∼ 3 keV. Finally, above m ∼ 300 eV, there is no resonance region in the Sun (298 eV is the highest photon effective mass in our solar model) and the HPs are always produced non-resonantly. These high mass HPs are produced mostly in the core, where most of the elements are fully ionised and therefore one does not see many atomic lines, with Iron being the obvious exception.

3D dynamic solar surface emission
So far we have assumed that the properties of the solar plasma, temperature, density and composition, are spherically symmetric. The solar interior is radiative and spherically symmetric to a good enough precision but, from r ∼ 0.7R Sun the outgoing energy transport becomes convective and inhomogeneities start to grow towards the solar surface. Energy is transported by hot upflows that expand progressively as they ascend, cool radiatively in the photosphere, and descend in a turbulent dense downflow. Hot diverging upflows bounded by cool lanes of downflows form the characteristic pattern of solar granules observed in the solar photosphere, see [104] for a review. The density and temperature profiles do not only depend on the solar radial coordinate. We have a full time-dependent 3D-problem.
Notwithstanding the complexity, the situation for the production of HPs has not changed much. HPs can be produced resonantly around a region where m 2 γ (ω, x) = m 2 and for frequencies in the visible and low mass HPs, this region is in very close to the solar surface. When we considered the Sun as spherically symmetric, the resonance region was a perfect spherical shell. Now we consider the inhomogeities due to surface granulation and we expect to have a textured surface which displaces slightly out of a perfect sphere with the hot upflows (because ionisation is larger) and in with the downflows. Solar granules have a characteristic extent of the order of 1 Mm, and evolve in times of the order of 10 minutes. These magnitudes are very large compared with the width of the resonance region and therefore it is reasonable to assume that our resonant flux equation (2.40) still holds locally. Since now the resonance region is not perpendicular to the radial vector we have to substitute the derivative |dm 2 γ /dr| by the gradient |∇m 2 γ | and average over the surface. We thus propose where the integral is over the surface (or surfaces) where m 2 γ (x) = m 2 . We have dropped the term involving e −τ * because resonances in the visible lie typically in the optically thick Sun.
We have made use of the 3D solar atmosphere model of [87] to estimate this 3D flux and compare it with the 1D idealisation shown before. The model covers a small 6Mm×6Mm patch of the surface, extends from -0.8 Mm to 0.5 Mm in depth and a period of 45 solar min, and thus gives a sensible description of the surface in length and time scales where granulation evolves. It was used to interpret the solar spectrum and determine atomic abundances from the shape and intensity of absorption lines. In general, the agreement between the fits and observations is excellent, which builds strong confidence in this model as a trustable description of the solar surface, at least in the line forming region. The model was used in [88] to assess the effect of HPs in the determination of solar abundances, which was found to be negligible. The model does not include magnetic fields and thus includes no sunspots or faculae. Since they are relatively scarce, we can neglect them.
A 2D time-slice of the situation is depicted in Fig. 17. We see the corrugated shape of the resonance locus for several HP masses (black lines) with a temperature background in colours. The function m 2 γ (ω, x) follows very closely the temperature because the ionisation fraction is most sensitive to it. HP masses below 10 −3 eV have resonance regions which are indistinguishably close in the picture to the region where m 2 γ (ω, x) = 0.  Figure 17. UP: 2D slice of the solar atmosphere model of [87] showing the temperature profile in eV as a function of depth and transverse extent along the solar surface (x). Superimposed are the resonance loci, m 2 γ (ω, x) = m 2 , where photon→HP conversions are resonant for m = 10 −4 , 10 −3 , 3.2 × 10 −3 , 4.5×10 −3 , 5.5×10 −3 eV from bottom to top. DOWN: For comparison we show the temperature profile of the spherically symmetric solar model of [89]. Note that the reference level to measure the depth is different in both cases.
We have estimated the integral (4.8) over a range of HP masses and energies and show our results in Fig. 18. Due to computational limitations, we have not integrated over the whole 3D model but just over one central 2D slice and then extrapolated the results to 3D. The procedure is as follows. For each of the 89  where X = 6 Mm and the brackets denote a time average However, in the real 3D calculation we have to correct for the corrugation in the coordinate y as well as the fact that the gradient will in general be larger than the 2D version. In these integrals, the number of points where the 2D-gradient is really low and thus the contribution to the integral is large is relatively scarce. These contributions are mostly suppressed by the third component of the gradient in the 3D case, but since they don't contribute much already we simply neglect them and approximate (on average) the inverse 3D-gradient with the inverse 2D gradient. All in all, we use (4.10) We remark that estimating the 3D flux by a x−average over radial profiles (i.e. 1D gradients) leads to very erroneous conclusions. The reason can be understood looking at Fig. 17. Around x ∼ 1 Mm there are regions where m 2 γ depends mostly on the depth x but not on d, the 1D gradient goes to zero in some cases and the flux would blow up. In the 2D calculation, the resonance width is much more controlled because there comparatively much less points where the 2D gradient goes to zero.
We show the time-averaged flux as a solid line and the 1 − σ contours due to time variability, which shall go to zero when a larger surface is integrated. The results are compared with the spherically symmetric estimate, shown as a dashed black line. Since the 3D simulation does not extend very deep in the Sun, we cannot cover the resonance region for the highest masses m > 5 × 10 −3 eV in this study and our results stop around that mass. At large depths d 0.5 Mm, the inhomogeneities are relatively small, variations shrink and the calculations seem to converge to the ∝ m 3 trend observed in the 10 −2 − 0.1 mass range derived before with the 1D averaged atmosphere. The last points of the 3D data lie typically a factor of ∼1.5 above the 1D calculations and we expect that the 3D and 1D lines converge before m ∼ 0.2 eV because resonances would happen deeper than 1 Mm where temperature and density fluctuations are small and 3D and 1D profiles should agree very well.
In general, the results agree well with our previous calculation in order of magnitude and trend but there are some O(1) differences. For HP masses in the meV, the flux is larger by a factor of a few and for smaller masses the flux is lower by a similar factor, although with some frequency dependence. Actually, there is a very simple interpretation of the differences. In general, 3D models of the solar surface have temperature profiles with are shallower in the solar interior and steeper in the outside when compared with the averaged 1D versions, see for instance Fig. 13 of [104] and compare our Fig. 17 up and down. This translates into the average profile for m 2 γ near the solar surface. Therefore, resonance regions which lie deeper in the Sun are larger and those lying outside are smaller than in the averaged 1D model used and the same translates into the HP fluxes. Since decreasing the mass and the frequency, the resonance region gets displaced towards the solar exterior, the HP fluxes are expected to decrease for low masses and low energies and increase for relatively high masses and high energies. This is precisely what we observe in Fig. 18. The 1D averaged profiles might be very useful for solar modelling but the formula for resonant HP emission (4.8) is highly non-linear in the temperature and density and thus it is clear that using them introduces new uncertainties. A precise determination of the solar HP flux below m ∼ 0.2 eV thus requires a full 3D model of the solar atmosphere.
Finally, note that the fluxes increase slightly also because, due to corrugation, the resonance region is larger. We estimate this effect to be at most 40%. The corrugation also influences the angular distribution of the flux, discussed in section 2.6. For HP masses below 0.02 eV, the resonance width ∆ψ * gets broadened on average to a value ∆ψ * ∼ 0.1Mm/R Earth ≡ 0.14 arcseconds. This is of the same order than the widening of ∆ψ * that an experiment would see if measuring in the whole visible range at the same time because of the energy dependence of r * , see Fig. 14.

Summary and conclusions
In this paper we have reviewed the phenomenon of photon↔HP oscillations in a plasma, focusing on a stellar hydrogen-dominated medium. We have shown how to build the polarisation tensor (index of refraction) of such a plasma and we have performed calculations for a spherically symmetric standard solar model [89] and a 3D model of the solar atmosphere [87]. The most complete atlas of transverse HP emission to date is presented in Fig. 16, and it is corrected at low energies and HP masses by the 3D calculations shown in Fig. 18. Tables with all the results are available and can be obtained by contacting with the author. We have studied in depth resonant photon↔HP oscillations in the conditions relevant for the solar plasma, including the effects of inhomogeneities. We identified the regions of the Sun where resonant conversions take place integrated the HP production from these regions. At low energies and HP masses, this contribution dominates the flux. For this task we were helped by two facts that we have also discussed: 1) that resonance regions are very thin with respect to temperature and density scale heights inside of the Sun and 2) the integral of the HP production across a resonance region is independent of the value of the mean-free-path of photons around it, i.e. it is the same for optically thick and thin resonances (as long as the mean free path is non-zero and the temperature is reasonably constant across its width). The HP flux from resonance regions only, in the 1D and 3D models is presented in Fig. 11 and Fir. 18. In general the flux grows towards lower energies, until overcome by the production threshold. For very low masses, where the threshold is not of concern, our calculations are limited by ou solar refraction model. Here, resonances happen in the solar atmosphere, where atomic transitions of metals and molecules might be relevant but were not included. In the ultraviolet range, the flux is smaller and has a non-trivial energy dependence, cf. Fig. 13, which we have discussed at length. In the visible, the energy dependence is very smooth and the effects of hydrogen lines are small, cf. Fig. 12. Differences between the 1-D and 3-D models have been found to be substantial, cf. Fig.18. Further refined work is necessary to improve the precision here, because our 3-D models did not extend much inside the Sun. Our calculations are uncertain by a factor of ∼ 2 below m ∼ 0.01 eV which decreases up to 20% around m ∼ eV.
We have found the visible flux very similar to previous estimates [69,70], which were however based on shaky grounds. The current calculation settles, to the authors knowledge, all the issues left previously undiscussed.
The computed HP flux at visible energies can be used to search for these particles with helioscopes in regions of parameter space where they might arise in string theories [38,42] and even constitute the dark matter of the universe [24]. Such a future helioscope should be much more powerful than present-day CAST, SUMICO or SHIPS to be more sensitive than the solar precision constraint [47] or the absence of ionisation events in XENON10 [50], both based on the much larger flux of longitudinally polarised HPs from the Sun. Due to the relatively modest flux of transversely polarised HPs, the magnitude of such an experiment is considerably big, although perhaps not unrealistic. The International Axion Observatory (IAXO) [105] is a proposed axion helioscope that seems to fall in the required category, large aperture and low background [106], although it mainly aims at the X-ray energies relevant for axions and deploys a huge superconducting toroid, which would not be necessary to search for solar hidden photons.

Acknowledgements
The author is gratefully indebted to Pat Scott and Martin Asplund for making available their 3D atmosphere model for the results of this paper. The input and inspiration from illuminating conversations with S. Troiksky, J. Jaeckel and G. Raffelt cannot be underestimated. Finally, it is a real pleasure to acknowledge the help of Davide Cadamuro during the early developments of this project. This work was supported by the the Alexander von Humboldt