Point-Particle Effective Field Theory III: Relativistic Fermions and the Dirac Equation

We formulate point-particle effective field theory (PPEFT) for relativistic spin-half fermions interacting with a massive, charged finite-sized source using a first-quantized effective field theory for the heavy compact object and a second-quantized language for the lighter fermion with which it interacts. This description shows how to determine the near-source boundary condition for the Dirac field in terms of the relevant physical properties of the source, and reduces to the standard choices in the limit of a point source. Using a first-quantized effective description is appropriate when the compact object is sufficiently heavy, and is simpler than (though equivalent to) the effective theory that treats the compact source in a second-quantized way. As an application we use the PPEFT to parameterize the leading energy shift for the bound energy levels due to finite-sized source effects in a model-independent way, allowing these effects to be fit in precision measurements. Besides capturing finite-source-size effects, the PPEFT treatment also efficiently captures how other short-distance source interactions can shift bound-state energy levels, such as due to vacuum polarization (through the Uehling potential) or strong interactions for Coulomb bound states of hadrons, or any hypothetical new short-range forces sourced by nuclei.


Introduction
Nature is full of examples where small but massive compact objects (of linear size R) interact with and control the motions of lighter neighbours within a much larger surrounding domain (of size a R). Examples include nuclei and atoms, stars and solar systems as well as legions of others. For such systems familiar arguments (such as the multipole expansion) show that only a few features of the compact object are often relevant to understanding motions in their larger environment. This simplicity usually emerges once observables are expanded in powers of small ratios like R/a.
Effective field theories [1,2] are the natural language for exploiting this kind of simplicity, though these are usually only formulated in a second-quantized language with all species of particles represented by their respective quantum field. For instance two-body contact interactions between two species of particles in a fully second-quantized framework would be represented in terms of their respective fields by terms like g(ψ * ψ)(χ * χ) in an effective Lagrangian.
Our companion papers [3,4] explore how to formulate such effective theories using instead a firstquantized language for the heavy compact object, reserving the second-quantized language for the lighter particles with which it interacts. 1 In this mixed first-quantized/second-quantized (one-two) language, if the heavy (χ) particle is in a position eigenstate situated at x = 0 then the two-body contact interaction mentioned above instead has the form g(ψ * ψ) δ 3 (x). This kind of formulation would be appropriate when the mass of the compact object is sufficiently large. In such situations all information about the source enters observables through the boundary conditions that are implied for the light fields at the position of the heavy compact object; boundary conditions that are completely determined by the source's first-quantized effective action.
This type of one-two formulation can have several advantages. One of these is the more direct connection it provides to the study of particle motion within a central (e.g. Coulomb or gravitational) potential, for which many useful tools are known (particularly for bound states). In this they are complementary to a fully second-quantized (two-two) formulation, such as for NRQED or NRQCD [7,8], in which induced quantities -like the nuclear Coulomb potential or solar gravitational fieldarise as a resummation of a particular class of interactions that dominate in some limits. By contrast, in the mixed one-two framework such classical fields are included into the zeroth order description about which one perturbs.
Furthermore, relating the near-source boundary conditions to the source action takes the guesswork out of small-r boundary conditions, and shows in particular why linear 'Robin' boundary conditions are so generic at low energies (see also [9]). More generally, they show how to handle singular potentials (like V (r) ∝ r p with p ≤ −2) unambiguously, despite the generic absence in these cases [10] of smooth solutions at the origin.
The study in [3,4] considered both nonrelativistic and spinless relativistic particles orbiting the massive compact object, focussing in particular on unusual effects that arise if the compact source size, R, is small enough that relativistic kinematics is relevant for the matching problem to the interior physics of the source even for bound states whose total energy, ω, is nonrelativistic: m − ω m. This mixed relativistic/nonrelativistic regime occurs when mR v 1, where v ∼ Zα is the speed of the orbiting particle (whose mass is m). (Here we take the source charge to be Ze and α = e 2 /4π is the usual fine-structure constant.) In particular, for relativistic spinless particles an interesting regime was identified for which energy shifts of S-wave states due to the source's finite size scale as where the last factor is the S-wave Schrödinger-Coulomb wave-function at the origin |ψ(0)| 2 ∝ (mZα/n) 3 . Effects like this, scaling linearly with R, are unusual and so lead to the question of whether similar shifts occur for the spin-half electrons and muons that arise in conventional and muonic atoms. We here address this question by extending the discussion of [3,4] to spin-half systems, finding that although many of the features of the Klein-Gordon problem of [4] also carry over to the Dirac field studied here the scaling of (1.1) does not: the corresponding leading Dirac expression instead gives the standard result: At first sight this difference in scaling may seem surprising, since spin-dependent effects in orbital energies might be expected to be suppressed by v ∼ Zα leading one to expect Dirac and Klein-Gordon predictions to agree at leading order in Zα. Although this expectation is true for most observables, it proves not to be true when tracking finite-size effects because relativistic effects are not small at radii r ∼ R once R < ∼ Zα/m ∼ (Zα) 2 a B (where a B is the Bohr radius). Indeed the ratio of δω KG and δω D given above is of order Zα/mR, which is order unity for electrons (for which mR ∼ Zα even though both are separately small).
Along the way we show how to formulate the near-source boundary condition for fermions, and why these differ from those that arise for bosons. We identify how the couplings for two-body contact interactions run, even at the classical level, and how this running goes over to the running found in [3,4] in the non-relativistic limit. This running properly captures how effective theories can sometimes generate scattering lengths that are much larger than the size R of the underlying object, and corresponds to the first-quantized version of a similar discussion found in [8].
Another result from [3,4] carries over to fermions: the fixed point of the running is not at c s = c v = 0 for charged sources (for which Zα = 0). It turns out this nontrivial fixed point is precisely what is required in order for the fixed point to reproduce standard results for the Dirac equation in the presence of specific nuclear charge distributions. That is, when we compare the PPEFT approach to explicit solutions to the Dirac equation in the presence of a finite-size charge distribution, we find that matching produces contact interactions for the PPEFT that sit precisely at the infrared fixed point of the RG flow. This shows why energy-level shifts take on a particularly model-independent form (proportional to the charge-radius r 2 p = r 2 and higher moments [11,12] -see also [13]) in the special case where the nucleus is modelled as a specific charge distribution.
In what follows we specialize for simplicity to parity-preserving interactions and spinless compact central objects, and so strictly speaking the interactions we find suffice in themselves to describe finitesize effects in the He + ion or muonic states in even-even nuclei [12,[14][15][16][17]. The effects we find also apply to nuclei with spin (such as hydrogen) once the effective theory of the first-quantized source is supplemented by the extra interactions that a nuclear spin allows. (We intend to return to discuss spinning sources more fully in a later paper.) In Section 2 we set the stage by introducing the point-particle effective action in the context of Dirac fermions. In Section 3 we derive the boundary condition and the induced renormalization group running in the presence and absence of a Coulomb potential. This leads to the discussion of bound state energy shifts implicated by the boundary condition in Section 4. We discuss applications of PPEFT for fermions in Section 5 and conclude in Section 6. We discuss various technicalities in the appendix.

Action and field equations
To make things concrete we focus on describing a relativistic spin-half charged particle interacting with a small charged source. The system of interest consists of a 3+1 dimensional 'bulk' action coupled to a 0+1 dimensional 'point-particle' action representing the small source (e.g. the nucleus of an atom), where W indicates the integration is over the world-line, y µ (τ ), of the source. In the final equality L B and L p are both regarded as being functions of the bulk fields evaluated at an arbitrary spacetime point, x µ . L p is also a function of the 'brane-localized' position field, y µ (τ ).

Action and field equations
Taking the bulk dynamics to be QED with a fermion of charge −e, the bulk action becomes with D µ ψ = (∂ µ + ieA µ )ψ. This should be considered in the spirit of a Wilson action, and so in principle also includes an infinite series of subdominant local terms involving more powers of the fields and their derivatives (whose effects are not important in what follows). The point-particle action is similarly given by an expansion in these fields, for which (for a spinless, parity-preserving source) the leading parity-even terms are 2 where the over-dot denotes differentiation with respect to proper time, the coefficients c s , c v andh all have dimension length-squared and the ellipses indicate terms suppressed by more than two powers of length. Notice that terms involving more than two powers of ψ first arise suppressed by a coupling with dimension (length) 5 , and so are nominally subdominant to several terms involving only two powers of ψ but more derivatives than those written above. Specializing to the rest frame for a motionless source,ẏ µ (τ ) = δ µ 0 , with charge Q = Ze the bulk field equations become and This last equality trades the parameterh for the mean-square charge radius: r 2 p = r 2 of the source charge distribution usingh = 1 6 Ze r 2 p .

Bulk solutions
We seek solutions to the bulk equations with a motionless point charge situated at the origin. The Maxwell equation is straightforwardly solved for the given source by choosing A = 0 and electrostatic potential Here the first term is the usual homogeneous solution to the Poisson equation, normalized using the boundary condition at small radial distance, r = , corresponding to nonzero electric flux r= d 2 Ω n · E = Ze . (2.8) This boundary condition can be obtained by integrating the Maxwell equation over small Gaussian pillbox of vanishingly small radius r = . By contrast, the second term in (2.7) is the particular integral arising when solving div E = −∇ 2 A 0 = 1 6 Ze r 2 p ∇ 2 δ 3 (x). We wish to repeat the above arguments for the Dirac field, whose field equation is where the second line specializes to energy eigenstates, 3 ψ(t) = ψ e −iωt , and to gauge potentials of the form (2.7). The parameter c v tot denotes the total localized combination This implies ψ L and ψ R are related by Outside the source these equations become ( / D + m)ψ = 0 which (see Appendix B for a summary in the present conventions) for rotationally and parity invariant situations have solutions of the parityeven form 12) and parity-odd form . (2.13) Here U ± are the spinor harmonics that combine the particle's spin-half with orbital angular momenta = j ∓ 1 2 to give total angular momentum j = 1 2 , 3 2 , · · · . The functions f ± (r) and g ± (r) are found by explicitly solving the radial part of the Dirac equation in the presence of a potential A 0 (r). For a Coulomb potential with source charge Ze these radial equations are (see Appendix B for details) together with These have as their general solutions Here A ± and C ± are integration constants, M[a, b; z] = 1 + (a/b)z + · · · are the standard confluent hypergeometric functions, ω is the mode energy and κ and ζ are defined by 18) with κ real because of our focus on bound states: m > ω. The parity of the solution enters the above formulae only through the parameter K = ∓(j + 1 2 ) where (perversely) standard conventions match negative (positive) K to parity-even (parity-odd) states.

Fermionic boundary conditions and the point-particle action
The next step is to formalize the boundary conditions at the surface of a spherical Gaussian pillbox of radius r = , along the lines of what is done in (2.8) for the Maxwell field. We now show how these relate the constants c s and c v of the source action to the ratios g + /f + and f − /g − at r = . These boundary conditions are again obtained from the source action by integrating the equations of motion over the interior of the pillbox using the delta-function.

Source-bulk matching
That is, given the action whereĝ ab = g µν ∂ a x µ ∂ b x ν is the induced metric on the world-volume of the source and N = c s +ic v tot γ 0 is the Dirac matrix specified by the source action S p . Then the ψ equation of motion is so integrating over the small Gaussian pillbox, P , of radius centred on the source then gives (in the limit → 0 of vanishingly small pillbox) Here n µ is an outward-pointing unit normal to the pillbox so n µ dx µ = dr, and the integral of the mψ term vanishes as → 0. Our conventions on gamma-matrices in polar coordinates are given in Appendix A. For spherically symmetric configurations (in the limit where is much smaller than all other scales of interest) this implies the boundary condition Notice this boundary condition is trivially satisfied pretty much anywhere in the absence of a source, for a small enough pillbox. This is because no source means c s = c v = r p = 0 and ψ varies slowly enough to be approximately constant across the pillbox. In this case the integral over all directions for γ r on the surface of the pillbox gives zero trivially. The boundary condition on the Gaussian pillbox can be written as d 2 Ω B ψ( ) = 0 where (3.5) The dimensionless coefficientsĉ s = c s /(4π 2 ) andĉ v = c v tot /(4π 2 ) can be interpreted as the coefficients of a term in a 'boundary action' defined on the codimension-one surface of the Gaussian pillbox, The subscript on B is meant to emphasize that the constantsĉ a (and in general also the original couplings c i themselves) also must carry an implicit -dependence if physical quantities are to remain unchanged as is varied (more about which below). To see what these boundary conditions mean we write them out separately for ψ L and ψ R , leading to Notice that these can be found from one other by making the replacements ψ L ↔ ψ R together with (ĉ v ,ĉ s ) ↔ (−ĉ v , −ĉ s ). Acting on bulk solutions (2.12) and (2.13) and evaluating the angular integrations, these givê In what follows we determine c s (R) and c v (R) from several hypothetical UV completions for the structure of the source of size R, and then regard (3.8) as a boundary condition that selects the exterior solution appropriate for the source of interest. This emphasizes that it is only through boundary conditions like (3.8) that the physics of a specific source can influence the exterior solution, and so enter into physical observables.

RG evolution
The radius of the Gaussian pillbox, r = , where the boundary condition is not a physical scale, and so must drop out of predictions for observables (unlike the physical size, R, of the underlying source, say). In detail, this happens because any explicit -dependence arising in a calculations of an observable cancels an implicit -dependence buried within the 'bare' quantities c s and c v . Following the procedure of [3,4] (which in turn builds on [18]), we next determine what the -independence of observables implies for the -dependence of c s and c v .
First we establish what is needed to ensure physical quantities remain independent of . Boundary conditions like (3.8) affect observables by determining the ratio of the integration constants that arise when integrating the bulk field equations. For instance, writing the general solutions, (2.16) and (2.17), to the radial part of the Dirac field equation in the form 9) it is the two ratios C + /A + and C − /A − that are determined by a boundary condition like the specification of (g ± /f ± ) r= . Energy levels for states of either parity are determined by demanding the resulting value for the appropriate C/A be consistent with what is required for C/A by normalizability of the modes at infinity. Scattering amplitudes are similarly determined by C/A. It follows that physical predictions are -independent if c s ( ) and c v ( ) are chosen to ensure C/A is -independent for both parity choices. At some level (3.8) says it all. Rather than reading (3.8) as fixing f ± /g ± at a specific radius given known values of c s and c v we can instead read the equations as giving c s ( ) and c v tot ( ) for known functions f ± (r) and g ± (r). This means that the -dependence of the right-hand-side of (3.10) is simply given by the r-dependence of f ± (r) and g ± (r) using (3.9), with r = . Because C ± and A ± are r-independent the above conditions tell us what c s and c v tot must do to keep them also -independent. Our greatest interest is when is much smaller than the typical scale a of the external problem (such as the Bohr radius, for applications to atoms), and in this limit it suffices to use the leading small-r form of the solutions f ± and g ± when computing the -dependence of c s and c v tot . In this regime solutions are usually well described by power laws, with (3.9) reducing to For such solutions the choice of C ± /A ± controls the precise radius at which one of these solutions dominates the other one, and as a result the RG evolution of the couplings implied by (3.10) in this regime describes the cross-over between these two types of evolution.

Non-relativistic limit
We start by examining this running for parity-even states in the nonrelativistic limit, which corresponds to the evolution found in [3,4] using the Schrödinger equation. The radial equations for parity-even states are given by (2.14) which imply in the nonrelativistic limit (for which the energy and mass are approximately equal, ω m, and much larger than all other scales) it follows that g + f + /(2m) f + . Using this in the second of eqs. (2.14) and dropping subdominant terms gives the Schrödinger equation (in the presence of a Coulomb potential), with Schrödinger field ϕ(r) = f + (r). In this limit the Dirac spinor is approximately given by (3.12) so in the nonrelativistic limit the combination appearing in the source action is 13) where h = c s + c v is the coupling for the analogous effective Schrödinger contact interaction. Defining the quantity λ := 2m h tot = 2m h + 2π 3 Zα r 2 p , the nonrelativistic limit of the boundary condition (3.8) therefore is 3.14) in agreement with the boundary condition found for a Schrödinger field coupled to a source with Lagrangian density L p = −h ϕ * ϕ δ 3 (x) [3,4]. These references also show that restricting to s-wave ( = 0) configurations and using the small-r asymptotic form ϕ 1 (r) ∝ r and ϕ 2 (r) ∝ r − −1 implies that for small the evolution of h given in (3.14) satisfies the differential RG equation 15) in which the last equalities defineλ. The evolution ofλ evidently has two fixed points, atλ = ±1, and these respectively correspond to λ = 0 and λ = −4π . Comparing with (3.14) shows these forms for λ are equivalent to having ϕ(r) ∝ r 0 and ϕ(r) ∝ r −1 (i.e. r and r − −1 for = 0), showing the crossover described below (3.11).

Relativistic running when Zα = 0
A similar story relates the solutions f and g to solutions of the Klein-Gordon equation in the relativistic case, as is most easily seen in the absence of the Coulomb interaction (Zα = 0), as we now show.

Parity-even case
When Zα = 0 the first of eqs. (2.14) again gives g + as the derivative of f + : for a mode of energy ω. Using this in the second equation then shows f + satisfies the Klein-Gordon equation. This shows that the r-dependence of the ratio g + /f + is proportional to the ratio χ /χ for a Klein-Gordon field: But refs. [3,4] show (even for Zα = 0) that if we define the quantity for χ a general = 0 solution to the Klein-Gordon equation, thenλ := (λ/2π ) + 1 satisfies the RG equation 3.19) for small enough to use the small-r asymptotic solution for χ(r). Here ζ s := 1 − 4(Zα) 2 . As Zα → 0 it follows λ as defined in (3.18) again satisfies the RG equation (3.15). These considerations show that when Zα vanishes, if we define the quantity 20) for parity-even j = 1 2 states, thenλ + D := (λ + D /2π ) + 1 satisfies the same RG equation, eq. (3.15), as doesλ in the Klein-Gordon case. Notice that in the nonrelativistic limit we have λ + D → 2m(c s + c v ) in agreement with the Zα → 0 limit of (3.14).

Parity-odd case
A similar argument goes through for the parity-odd j = 1 2 states. Parity-odd states satisfy the radial equations (2.15) and so when Zα = 0 we have Repeating the arguments of the parity-odd case then shows that g − = χ satisfies the Klein-Gordon equation and so implies thatλ − D = (λ − D /2π ) + 1 satisfies (for small ) the same RG equation, (3.15) as do the parity-even and Klein-Gordon cases, provided we define

Flow patterns
The flow obtained by integrating (3.15) is given (for small enough that f and g are dominated by their near-source asymptotic forms) bŷ a flow that is shown in Fig. 1. In the first equality the integration constant is chosen using the initial condition λ ± D ( 0± ) = λ ± 0 , while in the second equality η ± = sign(|λ ± D | − 1) and the RG-invariant quantities ± are defined as the scales where theλ ± D approach zero (or diverge). Which of these one uses depends on whether the RG trajectory of interest has |λ ± D | greater than or smaller than 1. In either case ± is given explicitly by inverting the first equality of (3.23): As shown in detail in [3,4], the -independence of physical quantities implies they depend only on λ ± D ( ) and through RG-invariant quantities like ± .
-�� Figure 1. Plot of the RG flow ofλ ± D (as defined in the main text) vs ln / when Zα = 0. A representative of each of the two RG-invariant classes of flows is shown, and is chosen as the place whereλ = 0 orλ → ∞, depending on which class of flows is of interest.

For
± (though not so large as to invalidate the small-r expansion of the mode functions at r = ) the flow approaches the fixed point atλ / this implies c s and c v simply become independent of in this limit.
For small the flow emerges from the repulsive fixed point atλ . Consequently for small the couplings c s and c v evolve linearly with (as opposed to the naive quadratic behaviour expected on dimensional grounds): and The flow describes the transition between these two asymptotic states, and clearly no source coupling (c s = c v = 0) is an RG-invariant fixed point, and it is also RG-invariant to have c v = 0 while c s runs (corresponding to + = − ). As a concrete example, suppose matching to a UV completion were to give the predictions , (3.27) and so ± R requires (g s ± g v )R −4π/(m ± ω). Unlike for the nonrelativistic case there is always an ω for which this can be satisfied, but because ωR 1 this is only possible in the effective theory if g s ± g v is sufficiently large and has the right sign.
For general the running couplings arê (3.28) which has the right limits for both large and small . Consequently , (3.29) which shows how the flow for ± is towards constant c s and c v , asymptoting to limits renormalized relative their values at = R.

Relativistic running when Zα = 0
We repeat the analysis of Section 3.2.2 this time for the case Zα = 0 as is relevant to the Coulomb problem.

Parity even
The running in the parity even case is determined by equation (3.8). The small radius expansion of the mode functions f (2.16) and g (2.17) yields to leading order (3.30) The RG running can be found by calculating the derivative d(ĉ s +ĉ v )/d and after inverting (3.30) inserting 2ζ as a function ofĉ s +ĉ v : (3.31) Defining the quantityλ which has the solution .

Parity odd
Similarly to the parity even case we can write (3.8) aŝ Repeating the procedure of the previous subsection we then find the running to be Again, one can define the quantityλ which has the solution (3.39)

Fixed points
From the running equations (3.33) and (3.38), it is clear that there are fixed points whenλ + D = ±ζ, and whenλ − D = ±ζ. However, from the solutions (3.34) and (3.39), we see that the fixed points of λ ± D are coupled. The fixed point obtained in the limit → ∞ (which we call the IR fixed point) corresponds to be λ + D = +ζ and λ − D = −ζ, so that The UV fixed point is similarly defined as the limit → 0 and is given by λ + D = −ζ and λ − D = +ζ, so thatĉ For later purposes (when comparing to results for specific nuclear charge distributions) we remark that the IR fixed point implies the couplings c s evaluate at = R to The attentive reader may also be puzzled as to why the running for Zα → 0 does not coincide with the Zα = 0 running found earlier. The reason for this is the observation that the limits → 0 and Zα → 0 do not commute, due to the appearance of factors of 1/(1 − ζ) 1/(Zα) 2 within the hypergeometric functions that furnish the Dirac-Coulomb solutions. (Related to this, mode functions can asymptote to r p at small r where p ∝ (Zα) 2 , again displaying non-commuting small-r and Zα → 0 limits.) As discussed in later sections, this makes the evaluation of energy shifts for bound states for specific values for Zα and nuclear size R somewhat subtle, since care must be taken to work to a consistent order in small quantities.

Higher-order interactions
For some applications it is insufficient to work only to lowest order in the nuclear size, and so we pause here to classify some of the next-to-leading interactions according to their dimension: where the operators appearing in L n has engineering dimension (mass) n . In this notation L 0 + L 1 + L 3 represent the terms already written in (2.3), so we now enumerate the dimension-4 interactions. At this order the operators consistent with invariance under rotations, gauge transformations and C, P and T are E 2 , B 2 and 4 ψ γ 0 D 0 ψ. We therefore take where c t and the 'polarizabilities'h E andh B are new effective couplings having dimension (length) 3 . For instance the time derivative appearing in the last of these terms contains contributions to the Dirac equation that resemble a correction to c v by an amount δc v ∝ c t ω. For nonrelativistic bound states and for c v ∝ R 2 and c t ∝ R 3 such corrections look like mR 3 |φ(0)| 2 contributions to the energy shift, and so contribute to some of the subleading corrections discussed below. One can continue in this way to as high a dimension as one wishes. Notice that the first interaction to involve more than two Dirac fields -such as 'three-body' interactions, like c 3b (ψ ψ) (ψ ψ) δ 3 (x)arises once we consider effective couplings with dimension (length) 5 .
Effectively, we can parametrize the boundary condition as Any microscopic source physics can only influence parity-even physical observables through their contributions to the constantsĝ i , only a few of which are relevant to any given order in the small expansion parameters. This makes these parameters useful proxies for specific models of source physics, and their values are computed in Appendix B for several simple examples. Although quantities likê g 2 can be traded for parameters like c t and/or h E we do not pursue this connection explicitly here.

Bound-state energy shifts
With a view to computing nuclear-size effects on atomic energy levels we next turn to the implications source contact interactions have for the energy of states bound to the source. Our assumptions of rotation invariance in S p restricts us for simplicity to atoms with spherically symmetric nuclei. What we find also applies to nuclei with spin but must be supplemented by spin-dependent nuclear-size effects (such as nuclear-size effects for hyperfine splitting [11]).

Energy-shift calculations
Bound-state energies are computed by reconciling the implications for the integration constants, C ± /A ± , appearing in (3.9) (or in more detail (2.16) and (2.17)) as imposed by the small-r and large-r boundary conditions. At small r the relevant boundary conditions are (3.8), which we repeat here for conveniencê 1) and the implications of these for C ± /A ± -as found using (2.16) and (2.17) -must be consistent with normalizability at large r, which implies The black curve plots the right-hand side of eq. (4. vs Zαω/κ, with the zero of energy chosen to be the eigenvalue of the n = 2 and j = 1 2 states. Standard Dirac energy levels correspond to places where the plotted quantity vanishes, while finite-size effects of the source correspond to those energies for which (4.2) instead equals a specified nonzero (positive) value. The dashed curves show two approximations to (4.2) that provide useful analytic expressions for energy shifts. The blue (red) curve shows the single-pole (double-pole) approximation to (4.2), described in the main text. In order to better display the shape of these curves, for plotting purposes we use ζ = 0.9 (and so Zα ∼ 0.45) and for concreteness expand about the pole at n = 2.
In the absence of a source the Dirac energy eigenvalues are given by solutions to C ± /A ± = 0, which (4.2) shows is satisfied when ζ − Zαω/κ = −N with N = 0, 1, 2, · · · . This returns the standard Dirac energy eigenvalues where n = N + j + 1 2 = 1, 2, 3, · · · is the usual principal quantum number. In the presence of a finite-sized source we instead solve for ω by equating the right-hand side of (4.2) to the nonzero value of C/A obtained by fixing f /g using the boundary condition (3.8) at nonzero r = . In practice this is done in two steps: (i) computing the value of C/A implied from the microscopic physics of the source (as parametrized by S p , say); and (ii) solving (4.2) for ω as a function of nonzero C/A, given a known form for C/A. We next consider each of these steps in turn.

Solving for δω
Solving for δω = ω −ω N with given C/A requires no knowledge of source structure since the right-hand side of (4.2) is dictated purely by the known solutions to the Coulomb-Dirac equation. Although this is easily done numerically, there are also accurate analytic approximations that are very useful (particularly when tracking the dependence of the result on external parameters), which are summarized briefly here. Figure 2 plots the right-hand side of (4.2) against energy with the zero of energy chosen to be the Dirac energy eigenvalue for a point-like source corresponding to a particular whole number N . Also plotted are two approximate forms, corresponding to approximating Γ(−N + δz) (−) N /[N ! δz] in just the denominator (single-pole approximation) or in both the denominator and numerator (doublepole approximation). As the figure shows, because of the presence of a nearby pole in the numerator the first of these approximations turns out only to have a radius of convergence of order (1−ζ) ∼ (Zα) 2 and so is only of use for extremely small δz. Figure 3. A plot of the relative error made when computing δω/|ψ(0)| 2 for nonzero C/A using two analytic approximations single pole (black solid) and double-pole (blue dashed) to the right-hand side of eq. (4.2) as described in the main text. The plot's horizontal axis is mR, where m is the mass of the orbiting fermion and R is the size of the source. For plotting purposes we use Zα = 1/137 and compute the shifts to the parity-even j = 1 2 state with n = 2 assuming the source to be a shell of positive charge with radius R.

������ ����
The double-pole approximation turns out to be much better then the single-pole one (particularly given that the left-hand side, Q := −(C/A), of (4.2) turns out to be positive for small δz), and suffices for identifying the leading energy shift and its first subleading correction. This can be seen in Fig. 3, which compares the solution obtained for δω using these approximate formulae to numerical results. For the purposes of these comparisons the source is assumed to be a fixed charged shell of radius R, whose energy eigenvalues can be computed exactly, and the state whose energy is perturbed is taken to be a parity-even S state (similar results obtain for parity-odd states). The plots show that the error obtained when using the double-pole approximation is order (Zα) 2 out to mR < ∼ O(1), for reasons identified below when we seek to compute O(Zα) 2 terms.
Concretely, the double-pole approximates the right-hand side of (4.2) using the leading Laurent expansion near the poles of the Gamma functions,  (4.4). This allows δx (and hence also δω) to be solved for explicitly as where (because our later focus is on j = 1 2 ) we trade N for the principal quantum number, n = N + 1, and write 1 − ζ 1 2 (Zα) 2 . (The single-pole approximate differs from the above by taking the denominator to be unity, and only gives the leading contribution reliably in the limit mR 1, if R is the typical size of the source.) It is useful to extract the naive Coulomb wave-function at the origin from δω by writing where we write m 2 − ω 2 N = c n (Zαm/n) 2 and so which can be taken as unity for the leading and O(Zα) correction but not once order (Zα) 2 contributions are required. As we shall see, for O(Zα) 2 corrections (4.4) must also be revisited to include also subleading terms in δx.
To use the above formulae in practice we require an expression for how Q = −C/A depends on the properties of the source. If the UV completion were a specific classical distribution, ρ(r), of radius R then C/A would be fixed by demanding continuity of f /g between the exterior and interior solutions at r = R (examples of this are discussed in more detail below). In general, knowledge of C/A is equivalent to knowledge of f /g at some radius, since this is ultimately the only way the physics of the source influences exterior phenomena. What is required then is an explicit expression for C ± /A ± as a function of f ± ( )/g ± ( ). In principle this is obtained by taking the ratio of expression (2.16) and (2.17) for the exterior solution (for each parity) and solving the resulting equations for C + /A + and C − /A − . This is efficient and easy to implement numerically and once this is done f ± /g ± at r = can be traded for constants in the source action through boundary conditions like (3.8).
Analytic expressions 5 for the required relation for C ± /A ± can also be found when is small enough to justify keeping only the leading small-r asymptotic form for the confluent hypergeometric functions in (2.16) and (2.17). Specializing to states with j = 1 2 -i.e. parity-even (S) states f + and g + and parity-odd (P ) states f − and g − -since these are the states most sensitive to finite-size effects of the source, we find the leading small-form where f ± /g ± is evaluated at r = and X is defined by X := (m − ω)/(m + ω). As we shall see, it is the factor of 2 in the exponent of (2κ ) 2ζ that is responsible for the main differences between this Dirac case and the Klein-Gordon problem studied in [4] (for which instead (2κ ) ζs appeared). This factor has its origins in the spin-orbit coupling that mixes two different orbital angular momenta into each state having fixed j.

Leading and first-subleading energy shifts
For detailed studies of the influence of nuclei on atomic energy levels one expands all contributions to bound state energies as a dual series in the small parameters (Zα) 2 and m Zα ∼ /a B , where a B = 1/(mZα) is the Bohr radius and R where R 1 fm is a typical nuclear size. In practice, comparison with experiments on atomic energy levels requires both the leading contribution and its subleading O(mRZα) correction, and for electronic atoms (Zα) 2 corrections are also required since for R of order a Fermi these are comparable in size to (mRZα) corrections. Our purpose in this section is to identify as generally as possible how these terms depend parametrically on the properties of the source.
Although (4.9) is sufficient for some applications, a more accurate approximation turns out to be required in order to track the leading subdominant coefficients in this kind of expansion. Increased accuracy is required for bound-state calculations because nominally independent variables like κ and X become specific powers of Zα once evaluated at the lowest-order bound-state energies ω = ω N . For instance, using (4.3) in the definitions implies and so higher powers of these compete with powers of Zα arising elsewhere (such as from the expansion of ζ). Extracting a particular order in Zα is further complicated by the appearance of factors of (1 − ζ) −1 ∝ (Zα) −2 in the expansion of the confluent hypergeometric functions M[a, 1 − 2ζ; ρ], due to the singularity of M[a, b; z] as b approaches a nonpositive integer. We next identify the leading and subleading O(mRZα) and O(Zα) 2 contributions to the energy shift. To do so we use the exact expressions, (2.16) and (2.17), for the general Dirac-Coulomb solution and solve for the integration constants Q = −(C/A) in terms of f /g evaluated at r = = R, finding where (as before) X := (m − ω)/(m + ω) and and with x = Zα ω/κ whilex = Zα m/κ. These are to be evaluated at the lowest-order solution, x = x N = N + ζ, where N = n − 1 and ζ 1 + 1 2 (Zα) 2 for j = 1 2 states, and we work only to subdominant order in mRZα and (Zα) 2 . Eq. (4.11) agrees with (4.9) at lowest order in ρ, for which M[a, b; ρ] = 1.
Since ρ = 2κ N R ∝ mRZα working to fixed order in Zα allows us to expand M in powers of ρ, but when doing so must be careful about factors of 1/(1 − ζ) ∝ (Zα) −2 appearing in the coefficients of the hypergeometric series. Such terms only arise when b of M[a, b; ρ] is a negative integer and so only are a factor in Q 11 and Q 21 . Since all powers of ρ involve the factor mR our guiding principle when expanding in ρ is to keep terms involving only a single subdominant power of Zα. This also allows us to neglect all subdominant powers of 1 − ζ ∝ (Zα) 2 in any ρ-dependent terms. Using ζ 1 − 1 2 (Zα) 2 and x x N = N + ζ N + 1 one finds and while and

Parity-odd leading energy shift
We next turn to parity-odd j = 1 2 P states (for which K = +1). In this case following the same steps reveals the leading contribution to be where the entire contribution of source physics is through (4.23) with (as before) X = (m − ω)/(m + ω).
Dropping all subdominant powers of Zα (and for consistency restricting the source contribution to ξ f f 1 (mRZα) gives the leading parity-odd energy shift As usual, this is smaller than the parity-even result because it is suppressed by the spin-orbit coupling required to link the P states to = 0 orbital angular momentum.

Parity-odd subleading O(mRZα) energy shift
Even though small, for some special cases (such as the charged shell described below) it happens that f 1 = 2 3 and so the leading contribution to parity-odd states vanishes. Such cases are dominated by the subleading contribution, for which where we can use ξ f =f 1 (mRZα) +f 2 (mRZα) 2 in the explicitly written factors, but stop at ξ

Subleading (Zα) 2 energy shifts
This section computes the subdominant O(Zα) 2 energy shifts for parity even and parity odd cases. Because factors of mR do not accompany the subleading powers of Zα it suffices to drop all nontrivial powers of ρ from the get-go and instead focus on the subdominant powers of (Zα) 2 . Because of this we can evaluate Q directly using (4.9), which is repeated here for convenience This is to be expanded to order (Zα) 2 , using ζ 1 − 1 2 (Zα) 2 and Using the corresponding terms in the source expansion ξ g ĝ 1 +ĝ 3 (Zα) 2 then gives Q + for parity-even (K = −1) states as (4.33) and it is also necessary to refine the double-pole approximation, by keeping subdominant terms in the Gamma-function expansion: where the harmonic numbers (see also The first term agrees with the leading result found earlier, and to these can be added the subleading (mRZα) corrections found in eq. (4.21) above. Some implications of these formulae are explored in the next sections.

Examples
As ever, the power in using an effective action to describe the short-distance properties of the source lies in its generality. That is, coefficients like c s , r p and c v can be used to describe the leading contributions due to any localized source physics, provided only that this physics arises over small enough scales, R, to make an expansion in powers of R/a useful (where a is a typical macroscopic scale -such as the Bohr radius of an exterior orbit). This ensures the model-independence of parametrizing physical quantities like energy shifts in terms of these parameters. This section emphasizes this point by indicating how several kinds of microscopic source physics contribute to effective couplings in the source action, S p , and how the above expressions reproduce familiar results in specific instances.

Explicit charge distributions
Perhaps the simplest example of microscopic source physics that can be parametrized by S p is the situation where the source is an explicit static charge distribution, ρ(x), rather than a point charge. Examples of this form are studied in the literature, with sensitivity to source structure often estimated by tracking how energy shifts alter as ρ(x) is varied through a plausible range of configurations [12,14,15,19,20].

Relations to moments
The leading terms in the source-dependent energy shift in this case have been calculated by perturbing the interior solution around the Coulomb problem and are known 6 to be given by [12] where µ is the reduced mass (so µ → m in the infinite-source-mass limit used here) and 2) and F N R and F REL are given in terms of various charge moments in [12]. This result is model-independent inasmuch as the expression for the coefficients of the series are universal functions of these moments, and so with the energy shift due to various charge distributions just differing in the values these distributions predict for the moments themselves. This is a more limited sense of 'model-independence' than we use here, since the model-independence of the predictions of the effective action apply not just to static charge distributions, but essentially to any kind of source physics that is sufficiently localized. (This model-independence of EFT methods for atomic measurements are emphasized within the 2nd-quantized framework in [17,21]. ) We verify in Appendix B that for a general static charge distribution, ρ(x), the quantityĝ 1 that dominates how source physics appears in g + /f + is related to the rms charge density, r 2 p = r 2 , by which implies that the leading energy shift given by (4.18) becomes (5.4) as required for consistency with (5.1). On the other hand, the boundary condition (3.8) shows how the parameterĝ 1 is also interchangeable with one combination of c s and c v tot through 6) i.e. the infrared fixed point found in (3.42). Note the difference from the Schrödinger running where we found that h = 0 is a fixed point that parametrizes a trivial boundary condition. The subdominant (mRZα) contribution also provides a relation betweenĝ 2 and the higher moment r 3 (2) . Comparing (4.21) with (5.1) and using (5.3) shows Although we do not have a general proof of this result, we can verify it for specific charge distributions. These higher terms can be related to higher-dimension interactions -such as those of (3.44) -in S p , using matching conditions similar to (5.5), although we do not pursue this here.

Specific charge distributions
The detailed calculations done for specific charge distributions [12,19,20] provide useful checks on the higher-order terms, since these must agree on the series coefficients for the specific charge distributions studied. To provide this check we compute the couplingsĝ i andf i for various charge distributions in Appendix B, and we here use these in the above expressions for h eff to verify agreement where overlap is possible.

Spherical charged shell
The simplest such example is that of a charged shell, for which which is convenient since the interior solution can be solved exactly in closed form. (We have checked that our numerical results for this case agree with those of [20].) For this distribution the rms charge radius is r 2 p = R 2 and r 3 (2) = 16R 3 /5.
For the parity-even state the boundary parameters appearing in g + (R)/f + (R) work out to bê (5.9) while for the parity-odd state the analogous parameters arê Using these values to compute the leading and subleading (mRZα) and (Zα) 2 energy shifts then gives for parity-even states. Notice the correct result for r 2 p and the cancellation of the n-dependence (and agreement with) the second moment r 3 (2) for this distribution. This expression also agrees well with numerical evaluation (as illustrated in Fig. 3).
In this case, becausef 1 = 2 3 , the leading parity-odd energy shift vanishes, leaving a result that is smaller than would naively be expected. The energy shifts predicted by the parametersf i in this case are Both of these results also depend on n in the way indicated by numerical evaluation.

Uniform spherical distribution
A second go-to example is the case of uniform charge distribution, although in this case the interior solution cannot be computed in closed form. We have verified that our solutions agree in this case with the numerical results given in [20]. Analytic expressions for the series expansion for the energy shifts are also given in [19], and we have verified that our results agree with these (and with [12]) in this case. Evaluating the boundary condition g + (R)/f + (R) using the interior solutions returns the following valuesĝ (5.13) while the same calculation for the parity-odd states giveŝ 5.14) Used in the parity-even energy shift, these values return the leading and sub-leading results These agree with the coefficients given explicitly in [19]. The first two terms also agree with [12] since the rms radius is r 2 p = 3 5 R 2 for this distribution, while the second moment is r 3 (2) = 32 21 R 3 and so

Other applications
A point-particle effective action like S p can be used to parametrize any short-range source physics and so need not be limited to describing the effects of finite nuclear size. This section summarizes a few such examples.

Vacuum polarization
A standard contribution to atomic energy levels that also can be captured using S p is the contribution (or parts of the contribution) due to vacuum polarization. It is well-known that the effects of vacuum polarization on the field of a point charge, Ze, can be described by the Ueling potential [22], of the form (5.17) in which m is the mass of the particle circulating within the loop. Since the range of this interaction is of order R ∼ m −1 the electron and muon vacuum polarizations fall into the category of physical effects acting over much smaller distances than typical sizes of orbits in ordinary atoms. The same is true for the influence of the muonic vacuum polarization within muonic atoms (but because m e ∼ αm µ it is not true for the shifts on muonic atom energies due to electron vacuum polarization). Such a potential shifts the energy of atomic states with low angular momentum that sample the potential near the nucleus, by an amount that is proportional (in the Schrödinger limit) to the wavefunction at the origin: |ϕ(0)| 2 . Using the notation of earlier sections, the resulting energy shift has size 18) where m is the mass of the particle in the loop. Since the photon line of the vacuum polarization does not flip helicity the arguments of earlier sections imply that this leading energy shift is correctly captured (at order (Zα) 2 /m 2 ) in all low-energy observables through a contribution to the effective couplings in (

Strong interactions and anti-protonic atoms
When the particle orbiting a nucleus experiences the strong interaction (such as for a π − , K − orp) then it experiences a short-range (R ∼ m −1 π ) strong interaction with the nucleus in addition to the usual Coulomb interaction. These are often described in the literature in terms of explicit nuclear potentials, which though concrete introduce an element of model-dependence into the treatment.
For such situations a more model-independent approach is to use the contact interactions appearing in (2.3) to capture the effects of these strong interactions on energy shifts and nuclear scattering amplitudes. This has the advantage of using only the short range of the force to organize the calculation, and so allows the disentangling of effects that rely only on this from those that instead depend on the detailed form assumed for any hypothetical nuclear potential.
Ref. [4] shows how parametrizing these strong interactions in terms of the lowest-dimension contact interaction allows the derivation of a relation between the strong-interaction induced shifts in atomic energy levels and the scattering length for collisions with the nucleus, that reproduce the standard Deser formula [23] (derived using nuclear potentials in the 1950s).
The leading effects of the nuclear force on antiprotons in protonium [24,25] can similarly be captured through the contact interactions of (2.3), though for protonium the existence of a relatively quick annihilation channel reduces the practical utility of using measurements of the energy shifts to learn about the nuclear interaction. But because this annihilation can also be described in the effective point-source action through the addition of imaginary parts to the effective couplings c s and c v one use for S p in this case is to compute the dependence of the annihilation rate on the principal quantum number n for S and P states. Thinking of the annihilation rate as the imaginary part of the energy eigenvalue shows that this n-dependence should be the same as for the energy shifts found in earlier sections, and this indeed reproduces what is found when modelling annihilation using nuclear potentials [26].
The virtue of rederiving this result using S p is that the effective field theory shows why the result is robust, and not an artefact of model-dependent details.

Exotic interactions
A fairly obvious use for contact interactions in the point-particle action is to parametrize the effects of any hypothetical new forces acting between nuclei and electrons or muons, and in particular forces that differ in strength between these two (since these can be captured through species-dependent values for c s and c v , unlike for r p ). Indeed the observation that the existence of such short-range interactions could, in principle, explain the proton radius puzzle [16] has led to efforts to better understand their size [21] and to the proposal of exotic interactions of this type [27].

Summary
In this paper, we introduce the PPEFT of Dirac fermions using a first-quantized language for the heavy compact object and a second-quantized language for the lighter fermion with which it interacts. This formalism can be advantageous to the fully second-quantized framework in the limit of the compact object being much heavier than the light interacting particle, i.e. the heavy compact object can be regarded as being in a position eigenstate to first approximation.
This formalism was previously introduced for bosons [3,4] where it was found that energy shifts due to the finite size R of the source scale linearly in R which is unusual. This does not carry over to fermions, i.e. energy shifts scale as R 2 . The absence of such unusual energy shifts means that there is no additional term that could account for the proton-radius-puzzle.
Our PPEFT allows one to parametrize the currently measurable energy shifts and their leading corrections due to the finite size of the nucleus for hydrogen and muonic hydrogen. Other applications include parametrizing strong interactions between the orbiting particle and the nucleus and antiprotonic atoms as well as hypothetical new forces acting between nuclei and electrons or muons.
In general, energy shifts are found by comparing the ratio of integration constants, C/A, appearing in the mode expansions (2.16) and (2.17) for the radial solutions to the Dirac equation found in two ways. On one hand normalizability at large r implies C/A is given by (4.2), while on the other hand it is fixed by the boundary condition for the ratio of radial functions, f /g, evaluated at a small radius r = near where the small source intervenes. The expression for C/A given f /g| r= can be found either by working numerically with the exact mode functions, or analytically using (4.9) if is small enough that the mode functions are well-approximated by their small-r asymptotic forms.
The main contribution of the PPEFT construction given here is to express f /g at r = in terms of general effective couplings, such as c s and c v using the conditions given in eqs. (3.8). This leads to a low-energy expansion applicable to a generic source physics provided only that the size of the source is sufficiently small. In the explicit calculations presented here 'generic' is in practice restricted for simplicity to parity conserving and rotationally invariant sources, rather than considering different source models one at a time. Results for specific models of the source can then be found by evaluating c s and c v explicitly using the model, such as along the lines as was done in the text for specific charge distributions.
What sets the size of ? The above procedure works for boundary conditions provided at any small radius r = , provided that is much smaller than the applications of interest (such as the Bohr radius, for atomic examples) while also being larger than the actual size R of the source. The effective couplings -e.g. c s and c v -themselves also depend on in precisely the way required to ensure that physical quantities do not; an evolution computed for c s and c v explicitly in §3.2. Once c s and c v are specified by matching to a specific model at r = R, their size at larger r = is dictated by this evolution.
Finally, we give explicit formulae for energy shifts in the Dirac-Coulomb case as a double series in powers of mRZα and (Zα) 2 , given a similar expansion for the boundary conditions f /g of the form 1 Zα g + f + r=R =ĝ 1 +ĝ 2 (mRZα) +ĝ 3 (Zα) 2 + . . . , (6.1) and with 'plus' and 'minus' referring to positive and negative parity eigenstates. The parametersf i and g i can be determined directly from a particular model of the underlying source and can be traded for parameters in the effective Lagrangian parameters (like c s and c v , with higher orders also depending on their higher-dimensional counterparts). Given such a boundary condition we write the energy shift to electrostatic bound states in terms of an effective δ-function potential: where the effective coupling h ± eff is given order by order in (Zα) 2 and (mRZα) by: These expressions apply for generalf i andĝ i out to subdominant order mRZα and (Zα) 2 , and so suffice for modern comparisons with precision measurements. As such they provide a modelindependent description of source effects, allowing source effects to be efficiently parameterized when comparing modern measurements [28] with other precisions corrections, such as those of QED. Finally, we have verified explicitly that these expressions reproduce those in the literature when specialized to the special case where the source is modelled as an explicit charge distribution, and for comparison purposes give expressions for the leading values off i andĝ i for several simple models considered elsewhere. and I is the 2-by-2 unit matrix. The gamma matrices are defined to satisfy the Dirac algebra {γ µ , γ ν } = 2 η µν where η µν is the inverse Minkowski metric, given (in rectangular coordinates) by diag(− + ++). Similarly and ψ := ψ † β = iψ † γ 0 . The chirality projection matrices are As usual, the Pauli matrices satisfy (A.5) and so defining γ µν := 1 2 [γ µ , γ ν ] we have Consequently the spin parts of the boost and rotation generators are block-diagonal in this basis, since

A.1 Polar coordinates
Our conventions for spherical polar coordinates {r, θ, φ} are standard, with (as usual) x = r sin θ cos φ , y = r sin θ sin φ and z = r cos θ . (A.9) The differentials therefore satisfy .10) in terms of which the flat 3D metric is This last equality defines the normalized frame of basis 1-forms, e r , e θ and e φ , so that an orthonormal frame is given by e r = dr , e θ = r dθ and e φ = r sin θ dφ . (A.12) We implicitly work in a gauge with ∂ µ A µ = 0. For later use notice the inverse of (A.10) is The radial gamma matrices then are defined by .16) and (for completeness) with σ r := (σ x cos φ + σ y sin φ) sin θ + σ z cos θ = cos θ e −iφ sin θ e iφ sin θ − cos θ 19) and so (A.20) which also implies For future reference notice also that with the convention 0rθφ = + det e a µ = +1/(r 2 sin θ) the above imply γ 0r = −γ r0 = γ 0 γ r = σ r 0 0 −σ r and γ θφ γ 5 = ir 2 sin θ σ r 0 0 −σ r , (A.22) and so .23) Solutions to the Dirac equation, ( / D + m)ψ = 0 also solve (A.24) which is the Klein-Gordon equation supplemented by a spin term, whose explicit form is

A.2 Spinor harmonics
When solving the Dirac equation we define quantities having definite quantum numbers (j and j z ) for J and J z , leading to the following 2-component spinors .
We directly evaluate for the case of most interest: j = 1 2 . For this purpose we use the explicit forms x ± iy r , (A.27) in the definitions of the U ± 1 2 jz to find (A.29) which are also what is found explicitly by acting on U + 1 2 jz with the explicit matrix Similarly, acting with σ r on U − 1 2 jz gives For later purposes we also evaluate the spatial derivatives explicitly using σ · ∇ = σ k ∂ k = σ x ∂ x + σ y ∂ y + σ z ∂ z as well as σ · ∇f (r) = f (r) σ · ∇r = f (r) σ · r/r = f (r) σ r . This trivially gives in agreement with an algebraic evaluation.

B Dirac solutions
This appendix collects several exact and approximate solutions to the Dirac equation that are used in the main text.

B.1 Exterior (Coulomb) solutions
Bound states for the Dirac equation are found as usual by demanding that the boundary condition (normalizability) at infinity be compatible with the boundary condition at the origin.

Energy eigenvalues
If the boundary condition at the origin is the usual one (for which we discard the singular solution to the radial equation -see below) the energy eigenvalues are where j = 1 2 , 3 2 , · · · and the principal quantum number is defined by n = N + j + 1 2 = 1, 2, 3, · · · . We define ζ = 1 2 1 + 4j(j + 1) − 4(Zα) 2 or so ζ → 1 as Zα → 0 when j = 1 2 . This implies 3) The standard derivation shows that for N = 0 (that is, except for n = j + 1 2 ) each state with fixed n and j comes two-fold degeneracy corresponding to parity s = ±. The most famous example is N = 1 and j = 1 2 , which corresponds to n = 2 and j = 1 2 in which case the degeneracy is between the 2S 1/2 and 2P 1/2 states that get split by the Lamb shift. This two-fold degeneracy does not occur for N = 0, corresponding to the n = j + 1 2 states like 1S 1/2 (the ground state), 2P 3/2 , 3D 5/2 and so on. (Notice that here S, P and D do not strictly correspond to specifying but instead give the parity value s for the corresponding state.)

Parity Eigenstates
Normally atomic states are given as parity eigenstates, which involves combining ψ L and ψ R since the action of parity is We expect a unique solution for each choice of parity, j and j z quantum numbers, while the above just relates the radial functions for left-and right-handed fields to one another. The Dirac equation and so To identify the parity eigenstates we expand in terms of the spinor harmonics U + and U − of Appendix A and define the radial functions f (r) and g(r) using the following ansätze: where the superscript on ψ and subscripts on f and g are the parity eigenlabel p = ±. Using this in either of (B.6) gives the same conditions relating g and f . For the parity even states the relations are while for parity odd states these relations instead become (B.9) as used in the main text.

Coulomb-Dirac solutions
To solve the radial Dirac equations, (B.8) and (B.9), for general radius we introduce the two functions where ρ = 2κr and κ = √ m 2 − ω 2 . Some manipulation shows that these satisfy the following firstorder linear ODEs which hold for either sign of the parity quantum number. The parameter ζ is as defined in (B.2). The most general solutions to these equations are given as linear combinations of confluent hypergeometric functions M(a, b; ρ) = 1 + (a/b)ρ + · · · , thereby introducing a total of four integration constants. The Dirac equation imposes two relations between the four constants. Hence, we can express the solutions Q 1 and Q 2 as where K = ∓(j + 1 2 ) for states with parity ±1. A and C are the two remaining integration constants, and are chosen so that the function multiplying A is bounded as ρ → 0 while the function multiplying C diverges there.
The corresponding expressions for f and g are then given by Normalisation of the state for ρ → ∞ demands A and C must be related by 16) which follows from the the large-ρ form of the confluent hypergeometric functions M. When C = 0 this condition reproduces the energy eigenvalue given in (B.1). Alternative boundary conditions at r → 0 change the bound state energy levels (and any other physical implications) entirely by changing what they imply for A/C. As the above formulae attest, such alternative boundary conditions governing A/C can be imposed by demanding that the ratio f /g take a specific value at a particular radius r = . (For instance, for particles orbiting a known charge distribution that extends out to radius r = R, it is continuity of the internal with the external solution at r = R that imposes the required condition: where f out and g out are the Coulomb solutions described above, valid for r > R, and f in and g in are given by the solving the Dirac equation for the charge distribution for r ≤ R. The next sections provide several representative solutions for simple charge distributions.

B.2 Interior solutions for given charge distributions
This section collects several simple solutions appropriate to the interior for several kinds of charge distributions, and gives the approximate series solutions in the general case.

B.2.1 Charged-shell model
In this case consider an exactly solvable model of a charge distribution against which later results can be compared. The model assumes a charge distribution that makes up a spherical shell, with surface density σ. That is, where R is the radius of the shell, and the second equality assumes the total charge is Ze. The corresponding electromagnetic potential found by integrating Maxwell's equations then is The Dirac equation outside the shell is therefore sees only the Coulomb potential and so is the one whose solutions are given above. The solution inside the shell is essentially the free Dirac equation, though in the presence of a nonzero constant A 0 . That is, it is equivalent to (A.24), which now reads (B.20) where the spatial derivatives are D i = ∂ i while the time derivative (acting on an energy eigenstate) is This has as solutions the usual spherical Bessel functions A j (kr) + B y (kr) , (B.22) and B = 0 if we demand R be bounded at r = 0. Specializing to j = 1 2 the appropriate solutions are f + = A + j 0 (kr), f − = B − j 1 (kr), g + = B + j 1 (kr) and g − = A − j 0 (kr). Since f + and g − vanish at the origin it follows that g + and f − must vanish there and this is automatic because these only involve = 1. When evaluated at r = R then Finally, the Dirac equation says f + = (m + W )g + and This allows the boundary condition to be written where we use (sin x − x cos x)/(x sin x) = 1 3 x + 1 45 x 3 + 2 945 x 5 + · · · . To make contact with the series for in powers of (Zα) 2 and mRZα we evaluate at a bound-state energy and use (kR) 2 = (ω + m)R + Zα (ω − m)R + Zα Similarly, for the parity-odd case These imply g + /f + − 1 3 (W − m)R − 1 3 (Zα) in the parity-even case for both the nonrelativistic and relativistic limits, while f − /g − 2 3 mR in the nonrelativistic limit (mR Zα) while in the relativistic limit (for which Zα/R ω m) we instead find f − /g − 1 3 Zα.

Expansion coefficients
For comparison with the results for other charge distributions for use in the main text it is useful to quote the above results in terms of parametersĝ i andf i appearing in the expansion g + f + r=R = Zα ĝ 1 +ĝ 2 (mRZα) +ĝ 3 (Zα) 2 + · · · , (B.31) and m − ω m + ω f − g − r=R = 1 2n f 1 (mRZα) +f 2 (mRZα) 2 +f 3 (Zα) 2 + · · · . (B.32) With these definitions the above calculation shows that the charged shell predicts for the parity-even state we haveĝ

B.2.2 General charge distribution
Next evaluate the interior solution for a general distribution ρ(r) for r ≤ R by evaluating as a series in kR. This is generally sufficient since kR M RZα or Zα in the cases mR Zα and mR Zα. The goal will be to determine f /g at r = R as a function of the first few derivatives of ρ at r = 0.

Parity-odd states
For parity-odd states the functions f − and g − satisfy (B.9), which reads ∂ r g − = m − ω + eA 0 (r) f − and ∂ r f − + 2f − r = m + ω − eA 0 (r) g − , (B.66) which has the same form as did the parity-even case if we make the replacements f + ↔ g − , f − ↔ g + and ω − eA 0 ↔ −(ω − eA 0 ). This implies the solutions have the same form with g ± i ↔ f ∓ i as well as M + ↔ M − and ρ i ↔ −ρ i .

B.2.3 Uniform charge distribution
A special case of the previous section is the case of a constant charge distribution ρ = 3Ze 4πR 3 for r ≤ R , (B.68) and so represents the special case ρ 0 = 1 and ρ k = 0 for all k = 0. For this distribution the rms radius and the moment r 3 (2) are given explicitly by The complete electrostatic potential therefore is eA 0 (r) = Zα R −1 + 1 2 u 2 − 1 , (B.73)