Axion-modified photon propagator, Coulomb potential and Lamb-shift

A consistent renormalization of a quantum theory of axion-electrodynamics requires terms beyond the minimal coupling of two photons to a neutral pseudoscalar field. This procedure is used to determine the self-energy operators of the electromagnetic and the axion fields with an accuracy of second-order in the axion-diphoton coupling. The resulting polarization tensor is utilized for establishing the axion-modified Coulomb potential of a static pointlike charge. In connection, the plausible distortion of the Lamb-shift in hydrogenlike atoms is established and the scopes for searching axionlike particles in high-precision atomic spectroscopy and in experiments of Cavendish-type are investigated. Particularly, we show that these hypothetical degrees of freedom are ruled out as plausible candidates for explaining the proton radius anomaly in muonic hydrogen. A certain loophole remains, though, which is linked to the nonrenormalizable nature of axion-electrodynamics.


I. INTRODUCTION
That the path integral measure in Quantum Chromodynamics (QCD) is not invariant under an axial chiral U(1) A −transformation provides a clear evidence that this classical symmetry does not survive the quantization procedure. As a consequence of this anomaly, QCD should not be a Charge-Parity (CP)-preserving framework. However, with astonishing experimental accuracy, no CP-violation event is known within the theory of the strong interactions. This so-called "strong CP-problem" finds a consistent theoretical solution by postulating a global U(1) PQ -invariance in the Standard Model (SM), which compensates the CP-violating term via its spontaneous symmetry breaking [1]. While this mechanism seems to be the most simple and robust among other possible routes of explanation, it is accompanied by a new puzzle linked to the nonobservation of the associated Nambu-Goldstone boson, i.e. the QCD axion [2,3]. As a consequence, constraints resulting from this absence indicate a feeble interplay between this hypothetical particle and the well-established SM branch, rendering its detection a very challenging problem to overcome. Still, various experimental endeavours are currently oriented to detect this elusive degree of freedom or, more generally, an associated class of particle candidates sharing its main features, i.e. axionlike particles (ALPs). Some of them being central pieces in models which attempt to explain the dark matter abundance in our Universe [4][5][6][7][8], whereas others are remnant features of string compactifications [9][10][11][12].
The problematic associated with the ALPs detection demands both to exploit existing high-precision tech- * Electronic address: selym@tp1.uni-duesseldorf.de † Electronic address: Alina.Golub@uni-duesseldorf.de ‡ Electronic address: c.mueller@tp1.uni-duesseldorf.de niques and to develop new routes along which imprints of these hypothetical particles can be observed [13,14].
Descriptions of the most popular detection methods can be found in Refs. [15][16][17][18][19]. A vast majority of these searches relies on the axion-diphoton coupling encompassed within axion-electrodynamics [20]. Correspondingly, many photon-related experiments such as those searching for light shining through a wall [21][22][23][24][25][26][27][28] and the ones based on polarimetry detections [29][30][31][32][33] have turned out to be particularly powerful. Contrary to that, precision tests of the Coulomb's law via atomic spectroscopy and experiments of Cavendishtype have not been used so far in the search for ALPs, although they are known to constitute powerful probes for other well-motivated particle candidates [34][35][36][37]. In particular, these setups provide the best laboratory bounds on minicharged particles in the sub µeV mass range. Simultaneously, by investigating the role of ALPs in atomic spectra, one might elucidate whether the quantum vacuum of these hypothetical degrees of freedom may be the source for the large discrepancy > 5σ between the proton radius that follows from the Lamb-shift in muonic hydrogen versus the established value based on electron scattering and the Lamb-shift in ordinary hydrogen [38][39][40][41]. Various theoretical investigations have been put forward seeking for a satisfactory explanation for this anomaly [42][43][44][45][46][47], some of them including hypothetical scalar particles.
Against the background of these circumstances, it is relevant to derive modifications of the Coulomb potential due to quantum vacuum fluctuations of axionlike fields and to study their potential consequences. The former are encompassed in the corresponding vacuum polarization tensor whose calculation, however, is not a straightforward task as far as axion quantum electrodynamics (QED A ) is concerned. This is because it requires-first of all-a meaningful implementation of the corresponding perturbative expansion in this nonrenormalizable framework. In analogy to quantum gravity [48][49][50][51][52][53][54], the expan-sion in terms of the axion-diphoton coupling gives rise to an infinite number of divergencies that cannot be reabsorbed in the renormalization constants associated with the parameters and fields of the theory. Unless a similar amount of counterterms is added, this feature spoils the predictivity of the corresponding scattering matrix, preventing the construction of a consistent quantum theory of axion-electrodynamics. Hence, the perturbative renormalizability of this theory demands unavoidably the incorporation of higher dimensional operators. This, in turn, comes along with the presence of a large number of free parameters which have to be fixed from experimental data. It is, however, known that this formal aspect relaxes because many of these higher-dimensional terms are redundant, in that the utimate scattering matrix is not sensitive to their coupling constants [55][56][57][58]. As a matter of fact, all contributions of this nature can be formally eliminated while the effective Lagrangian acquires counterms which allow for the cancellation of the loop divergences. Besides, when working at a certain level of accuracy, only a finite number of counterterms is needed and the cancellations of the involved infinities can be carried out pretty much in the same way as in conventional renormalizable field theories.
In this paper, axion-electrodynamics is regarded as a Wilsonian effective theory parametrizing the leading order contribution of an ultraviolet completion linked to physics beyond the SM. Its quantization is used for determining the self-energy operator of the electromagnetic field with an accuracy of second-order in the axiondiphoton coupling. This result is utilized then for obtaining the modified Coulomb potential of a static pointlike charge. In connection, the plausible distortion of the Lamb-shift in hydrogenlike atoms is established. Particular attention is paid to a limitation caused by the nonrenormalizable feature of axion electrodynamics which prevents us from having a precise and clear picture of the axion physics at distances smaller than the natural cutoff imposed by the axion-diphoton coupling. In contrast to previous studies of the Lamb-shift involving minicharged particles and hidden photon fields, this property introduces an unknown uncertainty that cannot be determined, unless the ultraviolet completion of axion-electrodynamics is found. We argue that-up to this uncertainty-axionlike particles are ruled out as plausible candidates for explaining the proton radius anomaly in muonic hydrogen (H µ ). Parallely, spectroscopic results linked to a variety of transitions in hydrogen are exploited to probe the sensitivity of this precision technique in the search for ALPs. We show that, as a consequence of the mentioned feature, high-precision spectroscopy lacks of sufficient sensitivities as to improve the existing laboratory constraints on the parameter space of ALPs.
Our treatment is organized as follows. Firstly, in Secs. II A and II B, techniques known from effective field theories are exploited for establishing the vacuum polarization tensor within an accuracy of the second order in the axion-diphoton coupling. There we show that, de-spite the electrically neutral nature of ALPs, the polarization tensor closely resembles the one obtained in QED. The similarity is stressed even further in Secs. III A and III B, where various asymptotes of the polarization tensor are established and the general expression for the axion-Coulomb potential is determined. The latter outcome is presented in such a way that a direct comparison with the Uehling potential can be carried out. Also in Sec. III B, we derive the corresponding modification to the Lambshift and emphasize the problematic introduced by both the effective scenario and the use of atomic s−states. In Sec. III C we provide some estimates and discuss the advantages and disadvantages of testing ALPs via excited states in H µ . Finally, in Sec. IV we expose our conclusions, whereas in the appendix a brief discussion on the sensitivity levels associated with precision tests of the axion-modified Coulomb law via experiments of Cavendish-type is presented.

II. THE MODIFIED PHOTON PROPAGATOR IN AXION-ELECTRODYNAMICS
A. Effective field theory approach Axion-electrodynamics relies on an effective action characterized by a natural ultraviolet scale Λ UV at which the U(1) PQ −symmetry is broken spontaneously. It combines the standard Maxwell Lagrangian, the free Lagrangian density of the pseudoscalar fieldφ(x) and an interaction term coupling two photons and an axion. Explicitly, Here, f µν = ∂ µāν − ∂ νāµ stands for the electromagnetic field tensor, whereas its dual readsf µν = 1 2 µναβ f αβ with 0123 = 1. Hereafter, we use a metric with signature diag(g µν ) = (1, −1, −1, −1), and a unit system in which the speed of light, the Planck constant and the vacuum permitivity are set to unity, c = = 0 = 1. As the axion-diphoton couplingḡ = Λ −1 UV has an inverse energy dimension, this effective theory belongs to the class of perturbatively nonrenormalizable frameworks. In the following we will suppose that Λ UV is a very large parameter in order to extend integrals over the momentum components to an infinite-volume Fourier space. Whenever no problem of convergence arises, the integrals over the spacetime coordinates will also be extended to the whole Minkowski space.
We want to use Eq. (1) to derive radiative corrections up to second-order in the axion-diphoton coupling ∼ḡ 2 . For this, we shall use well-established effective field theory techniques which have been applied extensively within the context of quantum gravity [48][49][50][51][52][53][54] and chiral perturbation theory [59][60][61][62] (see also Refs. [63][64][65] for their application in nonlinear QED). In connection, we will suppose that Sḡ characterizes-at energies substantially lower than Λ UV -the leading order contribution of its UV-completion which is Lorentz and gauge invariant. Hence, all higher-dimensional operators which are consistent with these fundamental symmetries should be included in Eq. (1), implying that an infinite number of counterterms is necessary to cancel out the divergences linked to one-particle irreducible Feynman diagrams [66]. However, to extract quantitative predictions from the radiative corrections in the second-order approach in the axion-diphoton coupling, it is sufficient to incorporate the next-to-leading order term of Sḡ. Combining the described method with a dimensional analysis, we find that the renormalization of the self-energy operators in QED A should be handled by two local operators of dimension 6: Various local operators sharing both the symmetry of the theory and the same dimensionality can be found. However, it can be easily verified that all of them reduce to those given in Eq.
(2) through integrations by parts. The Wilson parametersb a,φ determine the strength of the contributions above. They might be determined by a matching procedure provided the UV-completion of Sḡ is known. Since we ignore the precise form of the latter, they will be considered as arbitrary.
It is worth remarking that the action resulting from the combination of Eqs. (1) and (2) cannot be considered as an ordinary action containing higher-order derivatives. Within a quantum effective theory approach, higherdimensional operators-like those exhibited in Eq. (2)-are suppressed by higher powers ofḡ, so that their consequences at low energies relative to Λ UV =ḡ −1 are tiny when compared with the effects resulting directly from Sḡ [66,67]. To all effects, they must be treated as perturbations, otherwise a violation of the unitarity takes place due to the occurrence of Pauli-Villars ghosts. This problem can be noted straightforwardly by looking at the photon propagator of the electromagnetic theory that results from combining the Maxwell theory with the higherderivative operator written in the first line of Eq. (2). 1 Up to an unessential longitudinal factor this reads Here, a Pauli-Villars ghost is identified with the pole at p 2 =m 2 gh = (ḡb a ) −2 , the residue of which has an 1 The origin of this electrodynamics dates back to the work of P. Podolsky [68]. For further developments see Refs. [69][70][71][72][73][74][75][76][77].
overall sign different from the corresponding photonlike case [p 2 = 0]. Because of this, the norm of a corresponding quantum state excited from the vacuum would be negative, thereby violating unitarity. 2 The term written in the second line of Eq. (2) leads to a similar problem, but with a diffent ghost massm 2 s = (ḡb φ ) −2 . Symmetry arguments are therefore not enough to obtain a well-behaved quantum theory of the fields involved in S = Sḡ + Sḡ2 + . . .. However, if the scattering matrix linked to its UV-completion is unitary, it is natural to expect that the one associated with S, covering quantum processes at lower energies Λ UV , is unitary too. This idea justifies the restriction given above Eq. (3). We remark that these "ghost-providing" contributions are redundant operators which can be dropped from the effective Lagrangian without changing observables [55][56][57]. Later on, this outcome is used to remove them conveniently while the effective Lagrangian acquires counterterms which allow for the cancellation of the divergences associated with the loops that are calculated here [see below Eq. (7)]. Now, to carry out the renormalization program, the set of "bare" quantities (m,ḡ,m gh ,m s ,ā µ ,φ) should be replaced by the respective renormalized parameters (m, g, m gh , m s , a µ R , φ R ). In connection, each term in S = Sḡ + Sḡ2 + . . . has to be parametrized by a renormalization constant so that the action to be considered from now on is Here, T. D. Lee and G. C. Wick proposed a scenario in which the ghost is promoted to a degree of freedom [78,79]. They argued that heavy massive "ghost" modes should get a width, and thus decay into on-shell particles. This way, avoiding the existence of nonpositive norm states in the Hilbert space. The main presciptions to deal with this indefinite-metric quantum field theory were established in the cited references and in Ref. [80]. While this scenario lacks of a formal proof for the unitarity at an arbitrary loop order, there are no theoretical evidences that manifest the contrary. A study of our problem by using the Lee-Wick formalism is beyond the scope of the current work.
the renormalized axion field φ R (x) and its bare counter- Z φ is the corresponding wavefunction renormalization constant. Any other bare parameter relates to its respective renormalized quantity following multiplicative renormalizations according tō It is worth remarking that in Eq. (4) the shorthand where the connection betweenḡ and g has been used. Noteworthy, within our second-order approximation in the coupling constant, no modification on the renormalized axion-diphoton coupling g can be expected, i.e. Z g would not deviate from its classic tree-level value Z g = 1.
Hence, from now on, no distinction between the renormalized and physical coupling is needed. However, we emphasize that its bare counterpartḡ is still subjected to a renormalization due to the wave functions renormalization constants Z 3 and Z φ .
At the quantum level, the dynamical information of the system described by Eq. (4) is rooted within the Green's functions. They can be obtained from the generating functional where j(x) and j µ (x) denote the external currents associated with the axion φ R (x) and the gauge field a µ R (x). The contribution in the exponent of Eq. (7) which is proportional to the parameter 1 ζ guarantees a covariant quantization of a µ R (x). At this point it turns out to be convenient to bring the renormalized Lagrangian in S R [see Eq. (4)] to a canonical form in which the terms linked to the Pauli-Villars ghost are dropped. This can be achieved by performing the following local field redefinitions within the path integral [see Eq. (7)]: and by keeping the accuracy to the order g 2 . As these transformations are linear in φ R (x) and a µ R (x), the associated Jacobian leads to a decoupling between the corresponding Fadeev-Popov ghosts and the fundamental fields. The equivalence theorem [55,57] generalizes this fact by dictating that no change is induced on the scattering matrix through shifts of this nature; still they modify the initial parameters of the theory. Indeed, in our problem the redefinition of φ R (x) leads to a kinetic term of the form The additional factor contained in this expression can be reabsorbed in S ct [see Eq. (4)] by redefining the wavefunction renormalization and Z m → Z m to reabsorb terms arising when transforming S ct . Therefore, apart from this unobservable effect, the Pauli-Villars ghosts have no result other than to remove the divergences that might arise from the respective oneparticle irreducible graphs.

B. Renormalized photon propagator in a modified minimal subtraction scheme
We pursue our investigation by determining the modification to the photon propagator D αβ (x,x) due to quantum vacuum fluctuations of a pseudoscalar axion field φ R (x). The use of Eq. (7) allows us to express the modified photon propagator D αβ (x,x) as We then expand Eq. (7) up to the order g 2 and insert the resulting expression into the formula above. As a consequence, the corrected photon propagator reads Here Π µν (y,ỹ) encompasses the expression for the unrenormalized polarization tensor [see Fig. 1] as well as counterterms that allow for the cancellation of the divergences associated with this loop. Analytically, it reads where ∆ µν (x,x). Here, the external wavy lines represent amputated photon legs.
where the shorthand notationδ p1,p2 ≡ (2π) 4 δ 4 (p 1 − p 2 ) has been introduced and which diverges quadratically as |q| → ∞. The regularization of Π µν (p 1 , p 2 ) is then carried out by using a standard Feynman parametrization where C = e 1 2 (γ−1) /(4π) 1 /2 and γ = 0.5772... is the Euler-Mascheroni constant. In this context, µ denotes a dimensionful parameter, i.e. the substracting point that follows when rescaling the renormalized axion-diphoton coupling g → gµ /2 in D−dimensions so that its mass dimension −1 is kept. Also, when going from four to D−dimensions, the Wilsonian parameters [see Eqs. (2) and (6) Now, we integrate over q and Taylor expand the resulting expression in . As a consequence, Eq. (12) becomes with ∆(s) = m 2 s − p 2 s(1 − s). Manifestly, the term associated with the factor −1 is singular as → 0. In contrast to QED, such a divergence cannot be reabsorbed fully in the wavefunction renormalization constant of the electromagnetic field Z 3 by enforcing that the radiative correction should not alter the residue of the photon propagator at p 2 = 0 [66,81]. We solve this problem, by choosing the counterterms in the following form Notice that the ratio of scales mg acts like a dimensionless coupling constant. Thus, the one-loop renormalized polarization tensor in a modified minimal subtraction scheme (MS) scheme reads Noteworthy, this expression satisfies the transversality condition p 1µ Π µν MS (p 1 , p 2 ) = Π µν MS (p 1 , p 2 )p 2µ = 0. Some comments are in order. Firstly, when the QED action is extended with those terms belonging to QED A , the expression for Z 3 −1 found in Eq. (16) 16) and (17), the Fourier transform of Eq. (10) is, up to an unessential longitudinal term, This formula constitutes the starting point for further considerations. In the next section it will be used to establish the axion-modified Coulomb potential.

C. Axion self-energy operator, renormalized mass vs physical mass
Our aim in this section is to determine the axion selfenergy operator. Its associated Feynman diagram is depicted in Fig. 2. This object encloses the way in which the quantum vacuum fluctuations of the electromagnetic field correct the axion propagator. To show this analytically, we expand the generating functional for the Green's function [see Eq. (7)] up to first order in g 2 . Once this step has been carried out, the resulting expression is twice differentiated functionally with respect to the axion source j(x) leading to Here, Σ(y,ỹ) comprises the expression of the unrenormalized axion self-energy operator as well as some possible counterterms. Explicitly, Next, we Fourier transform Σ(y,ỹ) and regularize its divergent integral via dimensional regularization as made in Sec. II B. However, in contrast to the case treated there, the associated divergence at → 0 is fully reabsorbed here in a renormalization constant whereas Z m does not deviate from its classic tree-level value [Z m = 1]. This feature extends beyond the oneloop approximation because the axion-photon vertex prevents the proliferation of self-interacting terms for ALPs containing no derivatives [67]. Despite this, the bare massm still is subject to a finite renormalization due to the axion wavefunction renormalization constant [see Eq. (5)], which does not deviate from the classical value Z φ = 1. Keeping in mind all these details, we find that-in a MS scheme-the renormalized axion self-energy opera-tor is given by (15)]. As we could have anticipated, this expression is independent of the renormalized axion mass. It is, perhaps, worth stressing that-in a MS−scheme-the square of the physical mass m 2 phy is the value of p 2 for which the real part of the two-points irreducible function: 4 vanishes. Whenever the subtracting parameter satisfies µ m phy exp[−32π 2 /(g 2 m 2 phy )], m phy ≈ m holds. Hence, the expression for the polarization tensor [see Eq. (17)] as a function of m 2 phy would not differ from the one given in terms of the renormalized mass.
At this point we find it interesting to make a comparison between Eq. (22) and the polarization tensor given in Eq. (17). To this end, it is convenient to reexpress the latter as follows: In this formula, κ MS (p 2 ) represents the only nontrivial eigenvalue of the polarization tensor [82,83], which-in the limit under consideration-turns out to be smaller than Σ MS (p 2 ) by a factor −2/3. Let us finally remark that, in addition to the axiondiphoton interplay, axion self-coupling [84] as well as effective interactions with electron, proton and neutron might occur [16,[85][86][87]. In such a case further one-loop contributions to the axion self-energy operator might arise. electromagnetic four-potential where D αβ MS (p 1 , p 2 ) is given in Eq.
At this point it is worth emphazising that, while the expression for π MS (p p p 2 ) [see Eq. (17)] is finite, its dependence on the subtracting point µ introduces an arbitrariness. To remove it, we consider the expression of the electrostatic energy between two electrons in momentum space U(p p p). It can be established easily by taking the integrand above, with q → e R . After multiplying the resulting expression by e R , we find where the screened charge e 2 scr (p p p) = e 2 R [1 + π MS (p p p 2 )] has been defined.
As we still have freedom of performing finite renormalizations, we can demand that π MS (p p p 2 ) vanishes as |p p p| → 0. Since the corresponding length scale |x x x| → ∞, whereas the axion self-energy operator [see Eq. (22)] reduces to The results obtained so far are summarized in Fig. 3, which displays the behavior of π R (p 2 ) [left panel] and Σ R (p 2 ) [right panel] as a function of p 2 /m 2 . In both panels the respective real and imaginary parts are shown in green and red, manifesting by themselves the nonhermitian feature of the polarization tensor and the axion self-energy operator. To support this numerical evaluation from an analytic viewpoint, we first determine an exact expression for the imaginary parts. To this end, we restore the i0−prescription [m 2 → m 2 − i0] in Eq. (28) and apply the formula ln(−A − i0) = ln(|A|) − iπ with A > 0. Explicitly, where Θ(x) denotes the unit step function. We note that for an on-shell photon [p 2 = 0], the imaginary part of Π µν R vanishes, which implies-according to the optical theorem-that the emission of an ALP from a photon accompanied by the radiation of another photon is forbidden. This fact agrees with the outcome resulting from an analysis of the corresponding energy-momentum balance.
The expression for the imaginary part of the one-loop self-energy operator Im[Σ R (p 2 )] [second line in Eq. (30)] coincides with the result found previously in Ref. [88] through a direct application of the cutting rules. Its onshell evaluation [p 2 = m 2 ] should allow us to determine the total rate of the decay process φ → γγ via the relation Im[Σ R (m 2 )] = mΓ φ→γγ , provided the optical theorem is valid. With accuracy to first order in g 2 , this formula is indeed verified because an expression for Γ φ→γγ -relying on the corresponding S−matrix amplitude-can be inferred directly from the corresponding neutral pion decay rate [see for instance Eq. (19.119) in Ref. [89]]. This fact evidences that the unitarity is preserved, at least within the second order approximation in the axion-diphoton coupling g.
Further asymptotic expressions of Eq. (28) are elucidated. For m 2 p 2 , we find that π R (p 2 ) approaches to Conversely, for |p 2 | m 2 , its asymptotic behavior turns out to be dominated by the following function Notably, when p 2 > 0 is timelike, the expression above gets an imaginary contribution Im π R (p 2 m 2 ) ≈ −g 2 p 2 /(96π), which coincides with the leading order term of Eq. (30) when the condition p 2 m 2 is considered.

B. Electrostatic potential and modified Lamb-shift
The first contribution in Eq. (26) can be integrated straightforwardly, leading to the unperturbed Coulomb potential a C (x x x) = q/(4π|x x x|). The second one, on the other hand, will be computed by using π R (p p p 2 ) rather than π MS (p p p 2 ). With all these details we write a 0 (x x x) = a C (x x x) + δa(x x x), δa(x x x) = q d 3 p p p p 2 π R (p p p 2 )e −ip p p·x x x . For evaluating δa(x x x) explicitly, it is convenient to integrate by parts in Eq. (28) and use an equivalent representation of π R (p p p 2 ) instead: Observe that, for applying this formula in Eq. (33), p 0 must be set to zero. Taking this into account, the integral over the momentum can be carried out with relative ease. After developing the change of variable u = 1/(1 − s) 1 /2 , the axion-modified potential turns out to be We remark that spurious contributions containing Dirac deltas δ 3 (x) have been ignored since the theory is predictive only for distances |x x x| g. Although the integral involved in this formula can be calculated analytically, we will keep it as it stands. Mainly, because it will allow us to establish compact expressions for the energy shifts that atomic transitions undergo.
Asymptotic formulae for the modified potential can be extracted from Eq. (35) without much efforts. For instance, at distances larger than the Compton wavelength λ ∼ m −1 of the axion, i.e. for |x x x| λ, the region u ∼ 1 dominates in the integral involved in Eq. (35), and the axion-modified Coulomb potential approaches to However, at short distances [|x x x| λ], the main contribution to the integral in Eq. (35) results from the region 1 u (m|x x x|) −1 , and the integrand can be approached by its most slowly decreasing function in u, which is ∼ ue −m|x x x|u . Consequently, This expression is independent of the axion mass. Observe that the distance |x x x| must satisfy the condition |x x x| g/(4 √ 3π), otherwise our perturbative approach breaks down. Incidentally, the corresponding energy scale µ p ∼ |x x x| −1 ∼ 4 √ 3π/g coincides-up to a numerical factor of the order of one-with the Landau pole linked to QED A : µ L ≈ 4 √ 6π/g [for details see Ref. [67]]. The distortion of the Coulomb potential due to ALPs [see Eq. (35)] allows us to infer the induced modifications in the spectrum of a nonrelativistic hydrogenlike atom. Since so far no large deviations from the standard QED predictions have been observed, we will assume that these energy shifts are very small and, consequently, apply standard time-independent perturbation theory. When considering a first order approximation, the energy shift δε follows by averaging the correction to the electrostatic energy δU(x x x) = −e R δa(x x x) [see Eq. (33) and the Feynman diagram depicted in Fig. 4] over the 0th order wavefunctions |ψ n, ,j . Explicitly, We wish to exploit this formula to predict a plausible axion Lamb-shift for the 2s 1/2 −2p 1/2 transition in atomic hydrogen. For this case, Eq. (38) leads to where r ≡ |x x x| and R n stands for a radial hydrogen wave function. In particular, Here a B = (αm e ) −1 is the Bohr radius with α = 1/137 the fine structure constant and m e = 0.511 MeV the electron mass.
Observe that the integration in Eq. (39) covers the region [0, ∞). However, since QED A does not provide a precise information about the form of the axion-Coulomb potential for distances smaller than g, the integral over r must be splitted ∞ 0 dr . . . = g 0 dr . . . + ∞ g dr . . .. In the following we will assume that the contribution from the outer region [g, ∞) dominates over the inner region [0, g], which we ignore. 5 As we will see very shortly, the yet 5 Strictly speaking, in accordance with the treatment applied in Sec. II [read also below Eq. (35)], the splitting of the integral should be carried out at a certain point d fulfilling the condition d g. However, in order to avoid uncertainties stemming from this additional parameter, we set d = g. Our corresponding results should be considered as order-of-magnitude estimates, accordingly. undiscarded values for g turn out to be much smaller than any characteristic atomic scale. With all these details in mind we integrate over r and arrive at Let us study the asymptotes of this expression. We first consider the case in which a B λ. Under this condition, the term of the integrand [1 + ua B m] 4 is dominated by u 4 a 4 B m 4 . The integral resulting from this approximation can be computed exactly. After a Taylor expansion in mg 1, we find the compact expression δε ≈ αg 2 96π 2 a 3 B ln (gm) + γ + 11 12 .
As in Eq. (14), γ = 0.5772 . . . refers to the Euler-Mascheroni constant. In the opposite case a B λ, the integrand in Eq. (41) turns out to be a function that decreases monotonically with the growing of u. It is then justified to approach it through its most slowly decreasing part which is ∼ u 3 e −mgu /(1 + ua B m) 4 . As a result, the corresponding integral can be computed analytically by using (3.353.1) in Ref. [93]. In the limit of mg 1, it allows us to approach Notice that, the equation above is a good approximation whenever the condition a B g holds. Moreover, although Eqs. (41)-(43) apply for ordinary atomic hydrogen, they can be adapted conveniently for studying the same transition in other hydrogenlike atoms. When hydrogenlike ions with atomic number Z > 1 are considered, for instance, the correction to the Lamb-shift will be given by Eqs. (41)- (43), scaled by the factor Z and a B → a B /Z. If a muonic hydrogen atom is investigated instead, a replacement of the electron mass m e by the reduced mass of the system m r ≈ 186 m e would be required.

C. Precision spectroscopy in Hµ µ µ and the proton radius anomaly
Before continuing with the physics of virtual ALPs, we will estimate the contribution to the Lamb-shift by a meson whose interaction with the electromagnetic field resembles the one exhibited by ALPs [see Eq. (1)], i.e. the neutral pion π 0 . We should however emphasize that its effect should be understood as a consequence of the quantum vacuum fluctuations of its constituent quark fields. When thinking of the axion as π 0 , m → m π = 135 MeV and the coupling constant g → α/(πf π ) turns out to be determined by α and the pion decay constant f π ≈ 92 MeV [66]. Observe that the corresponding value of g ≈ 4.97 × 10 −3 fm is two orders of magnitude smaller than the proton radius r proton ≈ 0.876 fm [38]. The corresponding correction to the 2s 1/2 − 2p 1/2 Lamb-shift in hydrogen atoms is δε = −1.07 × 10 −12 meV. This value turns out to be five orders of magnitude smaller than the experimental uncertainty |δε 2σ | = 2 × 10 −7 meV, established at 2σ confidence level [35,37]. If the previous evaluation is carried out by considering a muonic hydrogen instead [m e → m r ], we find that the correction to the energy due to the neutral pion field is δε = −6.81 × 10 −6 meV. Since this is five orders of magnitude smaller than the existing discrepancy between the experimental measurement and the theoretical prediction δε = 0.31 meV [38,42,45], virtual neutral pions are excluded as possible explanation for the proton radius puzzle. Now, we wish to investigate whether the Lamb-shift induced by quantum vacuum fluctuations of axionlike fields might cure this anomaly. To this end, we will evaluate Eq. (41) considering a mass region in which reliable results can be extracted. Within a pure spectroscopy context, this occurs for ALP wavelengths smaller or of the order of the Bohr radius of H µ , i.e. λ a µ with a µ ≈ 285 fm, otherwise interactions of other nature must be included. Correspondingly, we can formally explore ALP masses fulfilling the condition 1 MeV m. However, in the range 1 MeV m 100 MeV the axion-diphoton coupling g has been constrained severely from various results, including those dealing with electron beam fixed-target setups [see compilation of bounds in Refs. [91,92]]. Conversely, the sensitivities in experiments where ALPs masses 100 MeV m 10 GeV are probed turn out to be much weaker [white sector in the right hand side of Fig. 5]. A recent investigation based on electron-positron colliders has constrained g to lie below g < 10 −2 GeV −1 [92].
A numerical assessment of the axion-modified Lambshift has been carried out by considering this yet undiscarded region. The outcome of this evaluation is summarized in the left panel of Fig. 6. Observe that the energy shift has been plotted in the form of log 10 (−δε[meV]). The highest value achieved for δε ∼ −10 −6 meV corresponds to g ∼ 10 −2 GeV −1 and m ∼ 10 −2 GeV. Toward higher axion masses m ∼ 10 GeV and lower axion-diphoton couplings g ∼ 10 −6 GeV −1 the correction to the Lamb-shift tends to decrease significantly [δε ∼ −10 −14 meV]. Both estimates coincide with the values resulting from Eq. (42). The smallness of these outcomes as compared with the aforementioned discrepancy rules out the corresponding virtual ALPs as candidates to explain the H µ anomaly. As we have anticipated above Eq. (41), the chosen values for 10 −7 fm g 10 −4 fm are much smaller than r proton ≈ 0.876 fm.
Clearly, the previous statements cannot be considered conclusive as our estimation undergoes theoretical uncertainties arising from both the internal limitation of QED A at short distances as well as the finite proton size. The latter being closely related to the fact that the sstates penetrate the nucleus deeply even for the chosen g. This last problematic can be relaxed if transitions between excited states with nonzero angular momentum are considered instead. Although the problem of their relatively short lifetimes constitutes a major issue for their experimental investigation, their measurements seem to be apriori a realiable way to have a cleaner picture of whether a certain ALP is the cause of the aforementioned discrepancy or not. Inspired by these arguments, we consider as an example the 3p1 /2 − 3d1 /2 transition. In this case, the required wave functions are An adequate replacement of the radial wavefunctions in Eq. (39) by those above allows us to determine the corresponding modification of the transition energy: In the limit a B λ, the expression above is well approached by δε ≈ −αg 2 /(4374π 2 m 2 a 5 B ). The expressions associated with H µ can be read off from the previous one by replacing a B → a B /(186). Taking as a reference the undiscarded region used previously, 100 MeV m 10 GeV with g < 10 −2 GeV −1 , Eq. (45) has been evaluated. The result of this assessment is shown in the right panel of Fig. 6. As in the 2s 1/2 − 2p 1/2 transition, the highest energy-shift linked to the 3p1 /2 − 3d1 /2 transition in H µ arises from the combination of g ∼ 10 −2 GeV −1 and m = 10 −2 GeV. In such a case δε ∼ −10 −11 meV, which is five orders of magnitude smaller than the outcome associated with the axion 2s 1/2 − 2p 1/2 Lamb-shift. It is worth emphazising that, while the uncertainty introduced by the use of spherically symmetric orbitals could be circumvented as described above, the one linked to the physics at distances shorter than g remains.

D. Sensitivity to ALPs in high-precision hydrogen spectroscopy
In this section we want to investigate whether the current sensitivity in atomic hydrogen can improve the existing bounds on the axion parameter space. To do this we will analyse the respective energy-shift in both 2s 1/2 − 2p 1/2 and 1s 1/2 − 2s 1/2 transitions. An expression for the latter δε = ∞ 0 dr r 2 δU(r) R 2 1s (r) − R 2 2s (r) can be easily determined by taking into account the formula for R 2s (r) in Eq. (40) and the radial part of the 1s1 /2 -state: R 1s (r) = 2 exp[−r/a B ]/a 3/2 B . Explicitly, While the right panel in Fig. 7 shows the energy shift for 2s 1/2 − 2p 1/2 [see Eq. (41)], the one in the left depicts the result associated with the 1s 1/2 − 2s 1/2 transition [see Eq. (46)]. Both evaluations have been carried out by considering the region of the coupling 10 −6 GeV −1 < g < 10 −2 GeV −1 . In contrast to H µ , reliable predictions from high-precision spectroscopy in ordinary hydrogen require to deal with ALPs masses m 10 −5 GeV, corresponding to wavelengths λ a B . The highest mass shown in both panels [m = 10 GeV] has been set in order to preserve the perturbative condition mg 1. When comparing the energy shifts resulting from each panel in Fig. 7 with the corresponding experimental uncertanities [|δε 2σ | = 2 × 10 −7 meV for 2s 1/2 − 2p 1/2 and |δε 1σ | = 1 × 10 −6 meV for 1s 1/2 − 2s 1/2 transition [37]] we conclude that, in order to improve the current bounds on the axion parameter space, an enhancement in sensitivity of at least five orders of magnitude is required. It is worth remarking that this sensitivity gap also manifests in other high-precision experiments searching for potential deviations of the Coulomb's law as those of Cavendish-type. For further details we refer the reader to Appendix A. This lack of sufficient sensitivity in the context of ALPs is significant when taking into account that these setups have allowed for constraining severely the parameter spaces of other weakly interacting sub-eV particles, including paraphotons and minicharged particles. However, we should emphazise that-in contrast  Fig. 6, the region for the axion mass m is wider, covering from 10 −5 GeV to 10 GeV. As previously, 10 −6 GeV −1 < g < 10 2 GeV −1 and the energy-shift is given in the form of log 10 (−δε[meV]).
to our investigation-these particle candidates have been treated within renormalizable frameworks and, thus, the bounds have been established on dimensionless coupling constants. Likewise, we have already indicated below Eq. (16) that in the axion theory the quantity playing the corresponding role combines two unknown parameters ∼ gm. Hence, the axion mass m suppresses the limits that can be inferred for g.

IV. CONCLUSION
Within the effective framework of axion quantum electrodynamics, terms beyond the minimal coupling of two photons to a neutral pseudoscalar field have been used to renormalize the polarization tensor and the axion selfenergy operator. The former outcome was used to establish the photon propagator distorted by the quantum vacuum fluctuations of axionlike fields, a piece essential for determining the modification of the Coulomb potential induced by both virtual photons and ALPs. This result allowed us to evaluate the way in which atomic spectra could change. Particular attention has been paid to the 2s1 /2 − 2p1 /2 transition in hydrogenlike atoms as it might constitute the most natural way of verifying our predictions experimentally. Likewise, this sort of axionmodified Lamb-shift has been considered in attempting to explain the proton radius anomaly in muonic hydrogen. By contrasting the experimental result with our theoretical prediction, it was found that-up to the uncertainties caused by the nature of the transition and the internal limitations of axion-electrodynamics-ALPs can be excluded as plausible candidates for solving the aforementioned problem.
Our investigation has revealed explicitly that neither atomic spectroscopy nor experiments of Cavendish-type allow us to infer bounds that improve the existing con-straints on the axion parameter space. This fact contrasts with analogous outcomes linked to scenarios containing minicharged particles and hidden photon fields, in which both precision techniques have turned out to be particularly valuable [34][35][36][37]. The loss of sensitivity within the axion context is conceptually rooted in the nonrenormalizable character of QED A and manifests-at the level of the modified Coulomb potential Eq. (35)through the dimensionless factor ∼ gm. This ratio of scales accomplishing somewhat a role similar to the coupling strengths of the photon-paraphoton mixing χ and the parameter in the minicharged particles scenario. To a certain extent the described problem justifies the existing demand for new laboratory-based routes looking for ALPs [13,14] by using strong electromagnetic fields [15][16][17][18], e.g., those offered by high-intensity lasers [94][95][96][97][98][99][100][101][102][103][104].
Let us finally remark that the expression for Π µν R (p 1 , p 2 ) [see Eq. (28)] constitutes an essential piece for a more general class of polarization tensors which result when external electromagnetic fields polarize the vacuum [83,94].
Today, these setups also provide the best laboratory bounds on mini-charged particles in the sub−µeV range [36]. Here, we want to estimate the sensitivity of this type of experiments in the context of ALPs. To this end, we consider a setup containing two concentric spheres: an outer charged conducting sphere-characterized by a radius b-and an uncharged conducting inner sphere with a radius a. Only if the electrostatic potential follows the r −1 law, the potential difference between the spheres vanishes and the cavity is free of electromagnetic field. However, deviations from the Coulomb potential like those induced by loop corrections [compare Eq. (35)] could lead to a nontrivial relative voltage dif- where A 0 (s) is an arbitrary potential in which the charge of the pointlike particle must be set to unity [106]. Now, to determine γ ab resulting from QED A we insert the axion-modified Coulomb potential a 0 (s, q = 1) [see Eq. (35)] into the expression above. Notice that, similarly to the case analyzed in section III B, the integral over s must be splitted into two parts: The integral that remains in this formula can be calculated analytically by using (3.351.4) in Ref. [93]. Since both b and a are macroscopic quantities, the conditions a, b, b − a g hold and we can approximate the expression above by Notice that Eq. (A3) is independent of the radius a of the inner sphere. When considering the limit mg 1 we obtain which does not depend on the axion mass m either.
With γ ab to our disposal, we can proceed to estimate the sensitivity of this setup in the search for ALPs.
For such a purpose, we use the benchmark parameters of the experiment performed by Plimpton and Lawton [a = 38 cm, b = 46 cm] in which a margin for γ ab exists, provided it lies below |γ ab | 3 × 10 −10 [107]. By making use of Eq. (A4) we find g 6.7 × 10 7 GeV −1 . (A5) We emphasize that-as a consequence of the perturbative condition [see below Eq. (A3)]-this result applies for m 15 eV. Besides, it is trustworthy for axion wavelengths smaller than the typical length scale of the spheres ∼ 0.1 m, i.e., for axion masses m 10 −6 eV. We note that the constraint in Eq. (A5) for the mass region 10 −6 eV m 15 eV has already been discarded by combining the experimental outcomes of collaborations such as PVLAS and OSCAR [see Fig. 5]. Hence, the sensitivity in this experiment of Cavendish-type is not high enough to improve the existing bounds on the axion parameter space.
It is worth remarking that a more accurate version of this kind of experiments has been carried out by using four concentric icosahedrons [108]. For obtaining first estimates, they may be treated approximately as four concentric spheres. In contrast to the setup of Plimpton and Lawton, here a very high voltage is applied between the outer two spheres with radii d = 127 cm and c = 94.7 cm. The voltage difference is measured between the two internal ones, with radii b = 94 cm and a = 60 cm, which are uncharged. This setup allows us to infer bounds for the ALPs parameters via the ratio between the voltage differences: .
(A6) We insert Eq. (35) into (A1) and evaluate the resulting formula in the various parameters contained in Eq. (A6). As a consequence, we end up with where δ ≡ c − d and terms of the order of ∼ g 4 m 2 have been disregarded. Notice that, in contrast to Eq. (A2), the integrand above lacks terms involving ∼ e −mgu . This could be anticipated because the numerator of γ abcd [see Eq. (A6)] does not contain a potential evaluated at the surface of the spheres [compare with γ ab given above Eq. (A1)]. Consequently, when the condition md = d/λ 1 is satisfied, the asymptotic expression for γ abcd becomes independent of the axion mass and quadratic in g: Since this formula applies for axion wavelengths larger than the typical length scale of the experiment, the outcomes resulting from it can be considered reliable as long as the interactions between ALPs and plausible fields/matter existing outside of the external icosahedron are negligible. Next, the aforementioned experiment achieves a precision |γ abcd | 2 × 10 −16 [36,108]. Combining this value with Eq. (A8) we constraint g to lie below g 9.8 × 10 7 GeV −1 for m 10 −7 eV.
Noteworthy, despite the improvement in the experimental accuracy, the resulting upper limit turns out to be comparable to the one found from the results of Plimpton and Lawton [see Eq. (A5)] and so, no improvement is found as compared with the existing constraints. The lack of sensitivity is understood here as a direct consequence of the quadratic dependence of γ abcd on g [see Eq. (A8)].