Annihilation diagram contribution to charmonium masses

In this work, we generate gauge configurations with dynamical charm quarks on anisotropic lattices. The mass shift of and charmonia due to the charm quark annihilation effect can be investigated directly in a manner of unitary theory. The distillation method is adopted to treat the charm quark annihilation diagrams at a very precise level. For charmonia, the charm quark annihilation effect barely changes the mass, but lifts the mass by approximately 3–4 MeV. For charmonia, this effect results in positive mass shifts of approximately 1 MeV for and , but decreases the mass by approximately 3 MeV. We did not obtain a reliable result for the mass shift of . In addition, we observed that the spin averaged mass of the spin-triplet charmonia is in good agreement with , as expected by the non-relativistic quark model and measured by experiments.


I. INTRODUCTION
The hyperfine splitting of charmonia, namely is a good quantity to test the precision of lattice QCD studies when charm quarks are involved. There have been many lattice efforts on from both the quenched approximation and full-QCD simulation with dynamical quarks. A recent full-QCD calculation gives MeV after the continuum extrapolation [1], which is consistent with the experimental value MeV [2]. The latest calculation carried out by the HPQCD collaboration finds MeV at the physical point after considering the quenched QED effects [3]. Obviously, this result, with a much smaller error, still deviates from the experimental value about MeV. Actually, most of the lattice QCD calculations of do not consider the contribution of charm quark annihilation diagrams (or disconnected diagrams) to the charmonium correlation functions, which may be the major cause of the systematic uncertainty. Although this kind of contribution is expected to be highly suppressed due to the Okubo-Zweig-Iizuka rule (OZI rule) [4][5][6], its contribution to can be sizable. To be specific, the contribution of disconnected diagrams to the correlation function is of order based on a naive power counting of the strong coupling constant , while it is of order for the case of . Especially, since is a flavor singlet pseudoscalar, its coupling to gluons can be enhanced due to the anomaly of QCD, which may shift the upward in the similar sense of the origin of mass [7]. Furthermore, given the lattice prediction of the pseudoscalar glueball mass GeV [8][9][10], which is close to the mass, the -glueball mixing may introduce an additional mass shift of . Of course the light hadronic intermediate states also contribute in the presence of the light dynamical quarks. Therefore, the above effects should be taken into account if one wants to obtain a more precise theoretical determination of .
There are also a few lattice studies on the annihilation diagram contribution to the masses of and [11,12]. It is found the mass of is affected little by this effect, but the mass can be sizably changed. Thus it is very possible that the annihilation diagram correction to comes mainly from the shift of . However, no quantitative results have been obtained due to the large statistical errors yet. A more sophisticated lattice study can be found in Ref. [13], where the disconnected part of the pseudoscalar correlation function with respect to the four-dimensional distance is complicatedly parameterized by intermediate states and quite a few excited states. The major results are that the disconnected diagrams result in an approximate +2 MeV shift of in the quenched approximation and an approximate +8 MeV shift in the full QCD case. However, for the lack of sea charm quarks, the mass shift of charmionium states stemming from the disconnected diagrams can only be derived indirectly and the complicated parameterization of may introduce theoretical uncertainties.

1S 1P
In this work, we investigate directly the contribution of charm sea quark loops to the masses of charmonium states. To do so, gauge configurations with charm sea quarks are necessary, such that the theory is in a unitary manner for charm quarks. As an exploratory study, we generate gauge ensembles on anisotropic lattices with degenerate sea quarks whose mass parameters are tuned to be close to the charm quark mass (the effect of light u, d and s sea quarks is ignored temporarily for simplicity). In the presence of charm sea quarks, the full correlation functions of charmonia (including the connected and disconnected parts) have well-defined spectral expressions, from which the masses can be derived directly.
The key task is to perform a precise calculation on the disconnected diagrams involved. Technically, we adopt the distillation method to treat the charm quark annihilation effects [14]. On the other hand, a large number of statistics is mandatory to get the precise signals of disconnected diagrams. Therefore, for the computational expense to be affordable, we generate large gauge ensembles on anisotropic lattices with the spatial lattice spacing being larger than the temporal lattice spacing . Based on this numerical prescription, we would like to investigate directly the effect of charm quark loops on the masses of and charmonium states.
As we will have a large number of statistics for charmonium correlation functions, we can test the 'Center of Gravity' relation for both 1P and 2P states and try to give an estimate of the mass of . In the quark model picture, the 'Center of Gravity' mass is expected to be equal to the mass of the spin singlet [15] for charmonium states. With a large number of statistics, it may be possible to test this relation precisely on the lattice.

1P 2P
This paper is organized as follows: In Sec. II we will give a brief introduction to the lattice setups for the calculation and the distillation method. Sec. III is devoted to the lattice derivation of mass shifts of and charmonium states owing to the disconnected diagrams. As a byproduct, In Sec. IV we check the agreement (or deviation) between the 'center of gravity' mass and the mass of the spin-singlet state for and charmonium from the connected part of charmonium correlation functions. The summary can be found in Sec. VI.

II. LATTICE SETUP AND THE DISTILLATION METHOD
A. Lattice setup The gauge configuration ensembles are generated on anisotropic lattices. The gluonic action is chosen to be the tadpole improved version [9,10,16,17] where and are the plaquette variable and the rectangular Wilson loop in the plane of the lattice, respectively, with referring to the spatial direction and referring to the temporal direction. The tadpole parameter is defined through the spatial plaquette , and is tuned self-consistently (The tadpole parameter of the temporal gauge links is set to be as usual). The parameter is the bare aspect ratio parameter with and being the lattice spacing in the spatial and temporal direction, respectively. For fermions, say, charm quarks in this work, the action is chosen to be the anisotropic version of tadpole improved tree-level clover action [10,17,18] proposed by the Hadron Spectrum Collaboration where and the dimensionless Wilson operator reads β m For the given and the bare quark mass parameter , the bare aspect ratio for fermions in Eq. (2) is tuned along with to give the physical aspect ratio that is derived from the dispersion relation of the pseudoscalar ( ) meson where and are the energy and the mass of the hadron in lattice units, and is the lattice momentum .
In practice, the procedure of the tuning of parameters is very complicated, since the parameters are highly correlated with each other. We omit the details here and just present the final ensemble parameters in Table 1.
As an exploratory study, we ignore the effect of light quarks and generate gauge configurations with flavors of degenerate charm sea quarks on an anisotropic lattice with the aspect ratio being set to be . The lattice spacing is determined to be fm through the static potential and . The bare charm quark mass parameter is set from the spin average of physical and masses. This lattice setup is based on two considerations: First, it is known that the signals of gluonic operators are very noisy and demand a fairly large statistics. Second, we will use the distillation method to tackle the quark annihilation diagrams (see below).

B. Distillation method
Quark-antiquark mesons can be accessed by quark bilinear operators on the lattice. Theoretically, these operators can couple with all the meson states and possible multi-hadron systems which have the same quantum number. If one is interested in the lowestlying states, the quark field can be spatially smeared as are color indices, refers to the Dirac spin component, and is the smearing function at timeslice . Hadronic operators constructed out of smeared quark field can dramatically reduce their coupling with much higher states. In order to obtain , we adopt the state-of-the-art approach, namely, the distillation method [14], which is briefly described as follows.
The distillation method starts with solving the eigenvalue problem of the gauge covariant second order spatial Laplacian operator on each timeslice of each gauge configuration where the spatial Laplacian is expressed explicitly as where are spatial indices and is the stout-link smeared gauge link. After that, for a given Dirac index , each eigenvector is used as the source vector to solve the linear equation array where the matrix is the Dirac matrix of fermions on the lattice and is the corresponding fermion propagator. This procedure is repeated for all the Dirac indices , the timeslices and the number of the eigenvectors . Finally, multiplying the eigenvectors to each of the solution vectors one gets for given and the matrix elements in the subspace expanded by the eigenvectors, which are usually called perambulators. Based on these perambulators, we can define a propagating matrix being the desired smearing function, where is a matrix with each column being one of the eigenvectors for and the row number referring to all the spatial points and color indices. Obviously the propagating matrix can be viewed as the propagator of the smeared quark field  Annihilation diagram contribution to charmonium masses Chin. Phys. C 46, 043102 (2022) in the background gauge field , namely, Note that if , then the completeness of implies that , such that is actually the all-to-all propagator of the original fermion field . One can define the norm of the smearing function to measure the degree to which the quark field is smeared. Due to the translational and the rotational invariance, can be averaged over all the coordinates that satisfying . It is easy to see that, if all the eigenvectors are involved in defining , then the completeness and the orthogonalization require . If the number of eigenvectors involved is truncated to be , then is smeared from a delta function to a function falling off rapidly with respect to . Figure 1 shows the profiles of in data points at and 100, where damps faster when becomes bigger, and the dashed curves are naive fits using Gaussian function forms in the interval . We take in this study. Now consider meson operators constructed out of the smeared quark field, J PC which have the same quantum number . The connected part of their correlation function can be expressed as where is the characteristic kernel that reflects completely the structure and the properties of operator . The last equation of Eq. (15) is the generic expression for mesonic two-point functions in the distillation method, where the perambulator is universal and the information of the meson operators are completely encoded in . The disconnected part of the correlation function can be similarly derived in terms of and as follows ] . (17) Finally, the full correlation function of and is 1S

1P
The Eq. (18) is the basic expression we apply to calculate the corresponding correlation functions of and charmonia.
S U (2) There are subtleties in considering the contribution from the annihilation diagrams to the charmonium masses since there are two degenerate sea quark flavors in this work. This means there is a global flavor symmetry, named as the isospin symmetry in principle. Thus the connected part of the correlation functions are actually if the flavor wave functions of are taken as . Similarly the disconnected part should be since gluons are isospin singlets and couple only with states. Therefore, the full correlation function of the states is 0 ++ A special attention should be paid to the channel. The quark loops in this channel have nonzero vacuum expectation values and should be subtracted. For this we adopt a 'time-shift' subtraction scheme [10] by redefining the correlation function as where the vacuum constant term cancels. In practice, we choose . The price we pay for this is a little loss in the signal-to-noise ratio.

TION IN DIFFERENT CHANNELS
Γ As the first step, we compare the relative magnitudes of the contribution of the annihilation diagram to the connected counterpart in each channel, which is measured by the ratio for a specific channel labeled by |r|/a s = 0 Γ where the subscript '11' refers to the correlation function of the operator in the operator set of the channel. The explicit operators are tabulated in Table 2. The at the two quark masses are shown in Fig. 2 and Fig. 3 for channels with and . It is seen that the magnitudes of are of order : The 's of , and channels are much larger than those of the other three channels. This is understandable since the charm quark loops in and channels are mediated by two gluons to the lowest order of QCD, while the charm quark loops in other three channels are mediated by at least three-gluon intermediate states. On the other hand, the disconnected part of can be also enhanced by the QCD anomaly. Especially, the of is two orders of magnitude smaller and turns out to be approximately independent of .

Γ
The contribution of the disconnected diagrams to charmonium masses can be estimated as follows. For convenience, we drop the superscript ' ' temporarily in the following discussion. Usually the connected part of the charmonium correlation functions and the full correlation functions can be parameterized as Since , the ratio of the disconnect part to the connected part can be expressed as Intuitively, each mass term in the correlation function is given by the propagator in the momentum space, namely, . When the disconnected diagrams are considered, the propagator acquires a self-energy correction , which is independent of the states. Thus the full propagator contributing to each term of can be expressed as with the relations Based on the OZI rule and from the observation in Fig. 2 and Fig. 3, it is safe to assume and . Thus we have where has been defined and is assumed to be varying very slowly with respect to . The latter is obviously justified by the observation that is usually of order or even smaller. If we assume is insensitive to in the energy range of interest, we can take the approximation and , which imply Thus we have where

t R(t) t
The last term in Eq. (30) will die out when increases. Therefore, if shows up a linear dependence on in a time region, the contribution of the disconnected dia- grams to the ground state mass, namely, , can be extracted by the slope of in this region.
(1, 2) 9] R(t) [9,14] t ∈ [4,9] 2σ t ∈ [10,14] χ 2 /do f [4,14] From Fig. 2 and Fig. 3, one can see that, for the , and channels, does show linear behaviors in different time regions. Therefore, in these time windows, we can drop the last term in Eq. (30) and fit using linear functions. The fit results are shown in these plots by curves with error bands, where the bands in darker colors illustrate the fit time window (the values of are in the lattice unit in the context). The fitted slopes and the 's on the two ensembles (labeled as E1 and E2) are listed in Table 3 along with the in the fitting time window given in the table. It should be noted that the large errors of the ratio function come from the disconnected diagram which becomes very noisy beyond . In the channel, shows good linear behavior in the time interval (very close to a constant) on both ensemble E1 and E2, so we attribute the behavior of in the interval to be mostly statistical fluctuations. Actually, the data points in this interval deviate from the linear behaviors determined from the interval by less than . If we include the data points in the interval to the fit, the central value does not change by much with a mildly larger . Actually on E2 ensemble, all the data points in the interval converge to a single line. Obviously, the of the is tiny and consistent with zero within the errors. This can be understood as the result of the strong suppression based on the Okubo-Zweig-Iizuka rule (OZI rule) [4][5][6], since the two charm quark loops are mediated by at least three intermediate gluons.
For and , the mass shifts due to the charm annihilation effect are small but non-zero. This is understandable since the contribution of two-gluon intermediate states is suppressed to some extent due to the Landau-Yang theorem [19,20]. It is surprising that this kind of mass shift of is relatively large and negative. The reason for this is not clear yet. Intuitively, lattice QCD studies predict that the lowest tensor glueball mass is approximately 2.2-2.4 GeV [8-10, 21, 22], which is lower than the mass and should result in positive mass shift χ c2 if and the tensor glueball can mix.

R(t) t G(t) G(t) G(t)
The situation of the scalar and the pseudoscalar channels is a little more complicated. It is seen in Fig. 2 and Fig. 3 that the ratio functions in these two channels increase rapidly when is large. This observation signals that there may be lighter states than the ground state charmonium contributing to the correlation function in the scalar and pseudoscalar channels. Actually, lattice QCD calculations predict that the masses of the lowest scalar and the pseudoscalar glueballs are 1.5 -1.7 GeV and 2.4 -2.6 GeV [8-10, 21, 22], respectively. When the disconnected diagrams are considered, these glueballs can contribute as intermediate states to the correlation function . Therefore, more mass terms should be added into the expression of (Eq. (24)) to account for the possible contribution from glueballs. Consequently, Eq. (25) for these two channels may be modified to where is a parameter resembling the mass difference between the ground state glueball and the corresponding charmonium state. Of course, one-glueball or multi-glueball states (if they exist) can also appear in other channels, but Fig. 2 and Fig. 3 show that their contributions have no sizable effects and can be ignored in practice. We tentatively use the functions form of Eq. (32) to fit of the scalar and the pseudoscalar channels, with in this equation being approximated by a linear function. The fit results are illustrated in Fig. 2 and Fig. 3 by curves with error bands, where one can see that the function form describes the data very well in large time ranges. This manifests the efficacy and the necessity of the second term in Eq. (32). The mass shift of the ground state pseudoscalar charmonium is determined to be +3.0(1) MeV for E1 ensemble and +3.1(2) MeV for E2 ensemble, respectively.

R(t) R(t)
is converted into the value in physical units using the temperal lattice spacing GeV. This result is consistent with the previous result in Ref. [13] and the result MeV derived from the glueball-charmonium mixing mechan-  Using the lattice spacing GeV, we get the mass shifts as follows: Since the fitted is very small, we do not quote its value here but only say that the charm annihilation effects in are negligible. The mass shifts from E1 and E2 ensembles are compatible with each other in most channels except for the channel, where the deviation of values on the two ensembles is larger than . The reason for this discrepancy is not clear yet and should be investigated in the future.

P-WAVE CHARMONIUM
Now that we have calculated the perambulators of charm quarks on the two gauge ensembles of large statistics, we would like to take a look at the charmonium spectrum. In this work, we will focus on the lowest two states in and channels. In the previous sections we have shown that the inclusion of the disconnected diagrams does not change the masses of charmonium states much, so in the following discussions we only consider the connected part of the correlation functions. For each quantum number, we will first build an operator set, then calculate the correlation matrix of the operators in this set. After that, we will solve the generalized eigenvalue problem (GEVP) to obtain the optimized operators that couple most to specific states. What follows are the technique details.

A. Optimized operators through variational method
As we addressed in Sec. II.B, the distillation method automatically provides the gauge covariant smearing function for quark fields on each timeslice. To be specific, the smeared quark field can be obtained by , where the duplicated ŷ c(x) spatial coordinate means the summation over the space volume. Thus, the quark bilinear operators for meson states in this study can be built in terms of the smeared charm quark field . We introduce the spatially extended operators for each channel as where is the multiplicity of with , and with being the specific combination of -matrix that gives the right quantum number of each of the and state (the explicit expressions of 's are tabulated in Table 2), and being the gauge transportation operator from to . Obviously, is gauge invariant. Thus one can get the kernel function Because the disconnected part is considerably smaller than the connected part and has little effect on present discussion, the correlation function can be approximated by their connected part as where is the connected part r |r|/a s = 0, 1, 2, 3, 4, 5, 6, 7 In practice, are chosen to be spatially on-axis displacements with . The differences between the different operators can be monitored through the effective masses of the corresponding correlation functions of the operators of these operators in the channel are shown in Fig. 4 where is defined. Therefore, the differences in show how different the operators in this operator set are, and also result in the different -behaviors of the effective masses of the corresponding correlation functions in the early time region. In each channel, we focus only on the lowest two states; therefore, we select operators with to compose an operator set, which are expected to couple to the interme- diate states more differently, as manifested by the effective masses of these correlation functions. Based on the operator set in each channel, we first carry out the variational method analysis to obtain the optimized operators that couple most to specific states. That is to say, we calculated the correlation matrix of each operator set and solve the generalized eigenvalue problem to get the eigenvector that corresponds to the -th largest eigenvalue . It is expected that the operator couples most with the -th lowest state. Thus we have the correlation function of the optimized operator , whose effective mass function can be defined similarly to Eq. (39).
functions for for all the channels are plotted in Fig. 5 (ensemble E1) and Fig. 6 (ensemble E2). and are clearly separated from each other, but still have mild time dependence in the early time range due to slight contamination from higher states. In order to obtain the mass values more precisely, we perform two-mass-term fits to In each of the six channels and , the mass plateaus corresponding to the lowest two states are shown. For each correlation function, a two-mass-term fit is performed with the second mass term being introduced to take into account the residual contamination of higher states. The darker colored bands illustrate the fit results using the corresponding time window.
Annihilation diagram contribution to charmonium masses Chin. Phys. C 46, 043102 (2022) 043102-9 where the second mass term is adopted to account for the higher state contamination. and are coefficients of the first mass term and the second mass term respectively. The fit results with the fitted parameters are also plotted in Fig. 5 and Fig. 6 as curves with error bands, where the extensions of the curves illustrate the fit windows. For all the fits, the values of the per degree of freedom ( ) are around one and exhibit the fit quality. The fitted masses of , charmonium states are tabulated in Table 4, in which the experimental values are also given for comparison.
The non-relativistic quark model expects that for the same principal quantum number , the multiplet with and has a 'center of gravity' mass , which is the spin average mass of this multiplet and should be equal to the mass of the spin singlet counterpart (the derivation of the relation can be found in the Appendix A). For example, for the multiplets multiplet, is defined as As shown in Table 4, the masses of the states are determined very precisely with the statistical error being less than 1 MeV, so we can check the relation of Eq. (44). The difference between the and the mass of the spin singlet is sometimes called the hyperfine splitting of charmonium states, which is denoted by in this study. On the two ensembles (E1 and E2) in this paper, we obtain as

1P 2P
In other words, the relation Eq. (44) is satisfied for states with a high precision. For states, we have

2P 2P
There is a substantial deviation from Eq. (44). We are not sure of the reason for this deviation yet. It is possible that charmonium states do have a non-zero hyperfine splitting. It is also possible that the masses of states are not determined precisely enough. This should be explored in future studies.
It is interesting to note that the experimental results support the relation of Eq. (44) to a very high precision [2]. For charmonium systems, the states have been well established. According to the PDG data, the 'center of gravity' mass of charmonia is MeV, which is almost the same as the mass MeV. For bottomium systems, the and states are below the threshold and have very small widths. They are approximately stable particles and have direct correspondence to the according states predicted by the non-relativistic quark model. The MeV and also reporduce the mass MeV and the mass MeV.
signed to be the candidates of the charmonia, which are (the old ), and , while there is no experimental evidence of yet. The resonance parameters of these three states are 's width is large, while those of the other two states are relatively small. Using Eq. (44), the mass of can be tentatively estimated to be V. SOME CAVEATS We would like to point out the possible caveats of our lattice setup. As we addressed in Sec. II, the tadpole improved tree level clover action proposed by the Hadron Spectrum Collaboration was adopted for charm quarks in this study [17,18]. Our calculation shows that the hyperfine splitting is approximately 60 MeV for the two ensembles E1 and E2. This is obviously very far from the experimental value MeV, and signals large discretization errors. The Hadron Spectrum Collaboration also found this large discrepancy in its calculation of the charmonium spectrum [24], where MeV was obtained. This discrepancy might be attributed to clover term of the fermion action, which results in the spin-orbital and spin-spin interactions between the charm quark and antiquark in the non-relativistic approximation. It is found [24] that a larger coefficient of the spatial clover term can give a larger . On the other hand, even though the masses of the and charmonium states are determined very precisely in this work (see Table 4), it is obvious that the fine splittings between the states are also smaller than those of the experimental values. In the σ ts F ts ∝ σ · E ∆ HFS non-relativistic approximation, these splittings are due to the spin-orbital and the tensor interactions [15] (also in Appendix A). The under-determined splittings in this work may imply that the coefficient of the temporal clover term of the fermion action also has large discretization uncertainties. It may also because we ignored light u,d,s quarks. Based on these observations, it seems that the fermion action we adopt in this work may not be a very good anisotropic version for charm quarks. It is noted that another version of the clover fermion action on the anisotropic lattices [25,26] can produce a larger [27][28][29]. Of course, whatever the lattice setup is, the continuum limit should be the same, but the smaller the discretization uncertainties, the smaller the systematic errors are after the continuum extrapolation.

VI. SUMMARY
In this work, we generate gauge ensembles with two flavor degenerate dynamical charm quarks on anisotropic lattices. This lattice setup enables us to investigate the contribution of charm quark annihilation diagrams to charmonium masses in a unitary theoretical framework for charm quarks. The distillation method is adopted to realize both the smearing scheme of quark fields and the disconnected diagrams in calculating meson correlation functions. With large statistics, the effects of disconnected diagrams can be derived with a high precision.
For a given quantum number , the mass shift of charmonium masses due to the disconnected diagrams can be derived from the ratio of the disconnected part of the charmonium correlation function to the connected part. It is found that the ratios are relatively large for the scalar and the pseudoscalar channel, and the time dependence of implies that there may be a contribution from states lighter than the corresponding lowest charmonium states to the correlation functions when the disconnected diagrams are considered. This is not surprising since glueballs can appear as intermediate states in this case, whose masses are predicted by lattice QCD to be lower than charmonium states. By considering the contribution of glueballs, the mass shift of the ground state nP Table 4. The masses of 1P and 2P states derived from ensemble E1 and E2. The masses are converted into the values in physical units. The 'center of gravity' masses are the spin averages of the spin-triplet charmonium masses. The experimental results are also presented for comparison.
The relation is expected by the non-relativistic quark model. We observe that this relation is satisfied to a very high precision level for charmonia, namely MeV and MeV at the two charm quark masses used in this work. For the charmonia, we observe a large deviation from this relation. There are two possible reasons for this deviation: this relation actually does not hold for states, or the masses of the states derived in this work have contaminations from higher states. This should be explored in future studies. On the other hand, for the two charm quark masses, the mass splittings between the and charmonium we obtained are smaller than those of the non-relativistic quark model predictions [30,31], but consistent with experimental results, given the assignment that , and are charmonium states.

ACKNOWLEDGMENT
Y. Chen also acknowledges the support of the CAS Center for Excellence in Particle Physics (CCEPP). The Chroma software system [32] and QUDA library [33,34]  In the quark model picture, the spectroscopy of charmonia can be studied by solving the non-relativistic Schrödinger equation of bound states with QCD-inspired inter-quark potentials along with the relativistic corrections. To order , where is the velocity of the (anti-)charm quark with mass within a bound state, the generalized Breit-Fermi Hamiltonian [15] is H =2m + is the spin-independent interaction term whose explicit form is irrelevant to the discussion here and is omitted, is the term for the spin-obital coupling, is the spin-spin interaction term, and is the tensor interaction part. If the potential can be split into the vector-like part and scalar-like part , the explicit expressions of , and are with S 12 ≡ 12 where are the spins of the charm quark and antiquark, and is their total spin. For a given state in terms of the eigenvalues of , the orbital momentum and the total angular momentum , the expectation value reads ⟨L · S⟩ = ] . (A5) As for the spin-spin interaction, if the is taken to be Coulomb-like, namely, , then is actually a contact interaction, whose expectation value vanishes for states, since their radial wave function at the origin is zero. For states ( and states), the spin-spin interaction results in the hyperfine splitting between the spin triplet and spin singlet states. The values of and for and states are listed in Table A1. It is easy to see that for a super-multiplet with , one has ⟨L · S⟩ avg ≡ 1 3(2L + 1)