text:
# 𝑉_{𝑐𝑠} from Pure Leptonic Decays of 𝐷_𝑠 with Radiative Corrections

## Abstract

The radiative corrections to the pure leptonic decay $`D_s\mathrm{}\nu _{\mathrm{}}`$ up-to one-loop order is presented. We find the virtual photon loop corrections to $`D_s\tau \nu _\tau `$ is negative and the corresponding branching ratio is larger than $`3.51\times 10^3`$. Considering the possible experimental resolutions, our prediction of the radiative decay $`D_s\tau \nu _\tau \gamma `$ is not so large as others, and the best radiative channel to determine the $`V_{cs}`$ or $`f_{D_s}`$ is $`D_s\mu \nu _\mu \gamma `$.

$`a`$, Department of Physics, FuJian Normal University, FuZhou 350007, China

$`b`$, Department of Physics, NanKai University, TianJin 300070, China

$`c`$, Institute of Theoretical Physics, Academia Sinica, P.O.Box 2735, BeiJing 100080,China

PACS numbers:13.20.Fc, 12.39.Jh, 13.40.Ks

The pure-leptonic decay $`D_s\mathrm{}\nu _{\mathrm{}}`$ can be used to determine the decay constant $`f_{D_s}`$ if the fundamental Cabibbo-Kobayashi-Maskawa matrix element $`V_{cs}`$ of Standard Model (SM) is known. Conversely if we know the value of decay constant $`f_{D_s}`$ from other method, these process also can be used to extract the matrix element $`V_{cs}`$. But there are the well known effect of helicity suppression we can see it by factor of $`m_{\mathrm{}}^2/m_{D_s}^2`$:

$$\mathrm{\Gamma }\left(D_s\mathrm{}\overline{\nu }_{\mathrm{}}\right)=\frac{G_F^2}{8\pi }\left|V_{cs}\right|^2f_{D_s}^2m_{D_s}^3\frac{m_{\mathrm{}}^{}{}_{}{}^{2}}{m_{D_s}^{}{}_{}{}^{2}}\left(1\frac{m_{\mathrm{}}^{}{}_{}{}^{2}}{m_{D_s}^{}{}_{}{}^{2}}\right)^2,$$
(1)
Of them only the process $`D_s\tau \nu _\tau `$ is special, it does not suffer so much from the helicity suppression, and its branching ratio may reach to $`4.5\%`$ in SM. However the produced $`\tau `$ will decay promptly and one more neutrino is generated in the cascade decay at least, thus it makes the decay channel difficult to be observed. For the channels $`D_se\nu _e`$ and $`D_s\mu \nu _\mu `$, besides the small branching ratios, there are only one detected finial state, the measurement of such channels are very difficult.

Fortunately, having an extra real photon emitted in the leptonic decays, the radiative pure leptonic decays can escape from the suppression, furthermore, as pointed out in Ref., with the extra photon to identify the decaying pseudoscalar meson $`D_s`$ in experiment from the backgrounds has advantages, since one more particle can be detected in the detector. Although the radiative corrections are suppressed by an extra electromagnetic coupling constant $`\alpha `$, it will not be suppressed by the helicity suppression. Therefore, the radiative decay may be comparable, even larger than the corresponding pure leptonic decays.

The radiative pure leptonic decays, theoretically, have infrared divergences and will be canceled with those from loop corrections of the pure leptonic decays. In all the existing calculation of radiative decays, this part is ignored, since they do not include the radiative corrections of the pure leptonic decays. In this paper, we are interested in considering the radiative decays and the pure leptonic decays with one-loop radiative corrections together. Since the process $`D_s\tau \nu _\tau `$ dose not suffer the helicity suppression and has a large branching ratio, the corresponding loop correction(virtual photon) to these process should has a considerable larger branching ratio, at least comparing with the radiative decay, and can not be ignored.

The contributions of the radiative decays are corresponding to the four diagrams in Fig.1. According to the constituent quark model which is formulated by Bethe-Salpeter (B.-S.) equation, the amplitude turns out to be the four terms $`M_i(i=1,2,3,4)`$:

$$M_1=Tr[\frac{d^4q}{\left(2\pi \right)^4}\chi (p,q)i\left(\frac{G_Fm_w^2}{\sqrt{2}}\right)^{\frac{1}{2}}\gamma _\mu (1\gamma _5)V_{cs}]\times $$

$$\frac{i\left(g^{\mu \nu }+\frac{p^\mu p^\nu }{m_w^2}\right)}{p^2m_w^2}ie[(p^{}+p)_\lambda g_{\nu \rho }+(kp^{})_\nu g_{\rho \lambda }+(pk)_\rho g_{\nu \lambda }]ϵ^\lambda \times $$

$$\frac{i\left(g^{\rho \sigma }+\frac{\left(pk\right)^\rho \left(pk\right)^\sigma }{m_w^2}\right)}{\left(pk\right)^2m_w^2}\overline{\mathrm{}}\frac{ig}{2\sqrt{2}}\gamma _\sigma \left(1\gamma _5\right)\nu _{\mathrm{}},$$
(2)
$$M_2=Tr\left[\frac{d^4q}{\left(2\pi \right)^4}\chi (p,q)i\left(\frac{G_Fm_w^2}{\sqrt{2}}\right)^{\frac{1}{2}}\gamma _\mu \left(1\gamma _5\right)V_{cs}\right]\frac{i\left(g^{\mu \nu }+\frac{p^\mu p^\nu }{m_w^2}\right)}{p^2m_w^2}$$

$$\times \overline{\mathrm{}}\left(ie\right)\overline{)}ϵ\frac{i}{\overline{)}k_{\mathrm{}}m_{\mathrm{}}}\frac{ig}{2\sqrt{2}}\gamma _\nu \left(1\gamma _5\right)\nu _{\mathrm{}},$$
(3)
$$M_3=Tr\left[\frac{d^4q}{\left(2\pi \right)^4}\chi (p,q)i\left(\frac{G_Fm_w^2}{\sqrt{2}}\right)^{\frac{1}{2}}\gamma _\mu \left(1\gamma _5\right)V_{cs}\frac{i}{\frac{m_s}{m_s+m_c}\overline{)}p+\overline{)}q\overline{)}km_s}\left(i\frac{e}{3}\overline{)}ϵ\right)\right]$$

$$\times \frac{i\left(g^{\mu \sigma }+\frac{\left(pk\right)^\mu \left(pk\right)^\sigma }{m_w^2}\right)}{\left(pk\right)^2m_w^2}\overline{\mathrm{}}\frac{ig}{2\sqrt{2}}\gamma _\sigma \left(1\gamma _5\right)\nu _{\mathrm{}},$$
(4)
$$M_4=Tr\left[\frac{d^4q}{\left(2\pi \right)^4}\chi (p,q)\left(i\frac{2e}{3}\overline{)}ϵ\right)\frac{i}{\left(\frac{m_c}{m_s+m_c}\overline{)}p\overline{)}q\overline{)}k\right)m_c}i\left(\frac{G_Fm_w^2}{\sqrt{2}}\right)^{\frac{1}{2}}\gamma _\mu \left(1\gamma _5\right)V_{cs}\right]$$

$$\times \frac{i\left(g^{\mu \sigma }+\frac{\left(pk\right)^\mu \left(pk\right)^\sigma }{m_w^2}\right)}{\left(pk\right)^2m_w^2}\overline{\mathrm{}}\frac{ig}{2\sqrt{2}}\gamma _\sigma \left(1\gamma _5\right)\nu _{\mathrm{}},$$
(5)
where $`\chi (p,q)`$ is Bethe-Salpeter wave function of the meson $`D_s`$; $`p`$ is the momentum of $`D_s`$; $`ϵ,k`$ are the polarization vector and momentum of the emitted photon.

As the $`D_s`$ meson is a nonrelativistic S-wave bound state in nature, higher order relativistic corrections may be small, being the first order approximation for a $`S`$-wave state, we ignore $`q`$ dependence in amplitude and in the $`D_s`$ wave-function and write the wave function of the meson $`D_s`$ as:

$$\frac{d^4q}{\left(2\pi \right)^4}\chi (p,q)=\frac{\gamma _5(/p+m)}{2\sqrt{m}}\psi \left(0\right).$$

Here $`\psi (0)`$ is the wave function at origin in the coordinate space, and by definitions it connects to the decay constant $`f_{D_s}`$:

$$\psi \left(0\right)=2\sqrt{m}f_{D_s},$$

where m is the mass of $`D_s`$ meson. Moreover we note that for convenience we take unitary gauge for weak boson to do the calculations throughout the paper.

There is infrared infinity when performing phase space integral about the square of matrix element at the soft photon limit. It is known that the infrared infinity can be cancelled completely by that of the loop corrections to the corresponding pure leptonic decay $`D_s\mathrm{}\nu `$.

If Feynman gauge for photon is taken, the amplitude of loop corrections corresponding to the box diagrams (a), (b) can be written as:

$$M_{\left(2\right)}\left(a\right)=\frac{2}{3}eA\frac{d^4l}{\left(2\pi \right)^4}\left[\frac{4i\epsilon ^{\alpha \mu \beta \nu }p_\alpha l_\beta 4\left(p_\mu l_\nu plg_{\mu \nu }+p_\nu l_\mu \right)+\frac{8m_c}{m_s+m_c}p_\mu p_\nu }{l^2\left(l^22plm_w^2\right)\left(l^2\frac{2m_c}{m_s+m_c}pl\right)\left[l^22l\left(pk_2\right)\right]}\right]$$

$$\times \overline{\mathrm{}}\left[2\left(pk_2\right)_\mu \gamma _\mu \overline{)}l\right]\left(\gamma _\nu \right)\left(1\gamma _5\right)\nu _{\mathrm{}},$$
(6)
$$M_{\left(2\right)}\left(b\right)=\frac{1}{3}eA\frac{d^4l}{\left(2\pi \right)^4}\left[\frac{4i\epsilon ^{\alpha \mu \beta \nu }p_\alpha l_\beta +4\left(p_\mu l_\nu plg_{\mu \nu }+p_\nu l_\mu \right)\frac{8m_s}{m_s+m_c}p_\mu p_\nu }{l^2\left(l^22plm_w^2\right)\left(l^2\frac{2m_s}{m_s+m_c}pl\right)\left[l^22l\left(pk_2\right)\right]}\right]$$

$$\times \overline{\mathrm{}}\left[2\left(pk_2\right)_\mu \gamma _\mu \overline{)}l\right]\left(\gamma _\nu \right)\left(1\gamma _5\right)\nu _{\mathrm{}},$$
(7)
where the $`l`$, $`k_2`$ denote the momenta of the loop and the neutrino respectively. These two terms also have infrared infinity when integrating out the loop momentum $`l`$.

After doing the on-mass-shell subtraction, the terms corresponding to vertex and self-energy diagrams (c), (d), (e), (f) can be written as:

$$M_{\left(2\right)}(c+d+e+f)=\frac{ieA}{4\pi ^2}\overline{\mathrm{}}\overline{)}p(1\gamma _5)\nu _{\mathrm{}}\times [ln\left(4\right)\frac{8}{9}+\frac{2}{9}\frac{m_sm_c}{m_s+m_c}ln\left(\frac{m_s}{m_c}\right)$$

$$+\left(\frac{2}{9}+\frac{8}{9}\frac{m_c}{m_s+m_c}\right)ln\left(\frac{m_s+m_c}{m_s}\right)+\left(\frac{8}{9}+\frac{8}{9}\frac{m_c}{m_s+m_c}\right)ln\left(\frac{m_s+m_c}{m_c}\right)$$

$$+\frac{2}{\epsilon _I}2\gamma +ln\left(\frac{4\pi \mu ^2}{m^2}\right)+ln\left(\frac{4\pi \mu ^2}{m_e^2}\right)].$$
(8)
where $`A`$ is:

$$A=\frac{f_{D_s}\left(\frac{G_Fm_w^2}{\sqrt{2}}\right)^{\frac{1}{2}}V_{cs}eg}{2\sqrt{2}}=f_{D_s}\left(\frac{G_Fm_w^2}{\sqrt{2}}\right)V_{cs}e.$$

The other loop diagrams(we do not show them) always have a further suppression factor $`m^2/m_{w}^{}{}_{}{}^{2}`$ to compare with the loop diagrams we considered, and there is no infrared infinity in these loop diagrams, we can ignored their contributions safely. Furthermore we should note that in our calculations throughout the paper, the dimensional regularization to regularize both infrared and ultraviolet divergences is adopted, while the on-mass-shell renormalization for the ultraviolet divergence is used.

Detail cancellation of infrared divergence is given in Ref. Here we simply show the results. The ‘whole’ leptonic decay branching ratios, i.e., the sum of the corresponding radiative decay branching ratios and the corresponding pure leptonic decay branching ratios with radiative corrections, and put them in Table (1). The reason we put the radiative decay and the pure leptonic decay with radiative corrections together is to make the branching ratios not to depend on the experimental resolution for a soft photon. For comparison, the branching ratios of each pure leptonic decay at tree level is also put in Table (1). The values for the parameters $`\alpha =1/132`$, $`|V_{cs}|=0.974`$, $`m_{D_s}=1.9686`$ GeV, $`m_s=0.5`$ GeV, $`m_c=1.7`$ GeV, $`f_{D_s}=0.24`$ GeV and the lifetime $`\tau (D_s)=0.469\times 10^{12}s`$.

Table (1) Branching Ratios of the ‘Whole’ and Tree Lever Leptonic Decays

|  | ‘whole’ | tree |
| --- | --- | --- |
| $`B_e\left(10^5\right)`$ | 2.56 | 0.0108 |
| $`B_\mu \left(10^3\right)`$ | 4.706 | 4.605 |
| $`B_\tau \left(10^2\right)`$ | 4.138 | 4.489 |

We can see that, the ‘whole’ decay branching ratios $`Br_e`$ and $`Br_\mu `$ are larger than the corresponding branching ratios of tree lever, while the ‘whole’ $`Br_\tau `$ is smaller than the tree lever one. It means the contributions of loop corrections are negative, the dominate contributions of first order corrections to the pure leptonic $`D_s`$ decays are radiative decays when the lepton is $`e`$ or $`\mu `$, and is loop corrections when the lepton is $`\tau `$. So, the loop contributions are important for the decays $`D_s\mu \nu _\mu `$ and $`D_s\tau \nu _\tau `$, especially for the later. Through Table (1), we obtained that the rediative decay has a branching ratio $`Br(D_s\mu \nu _\mu \gamma )>1.01\times 10^4`$ and the loop correction to $`D_s\tau \nu _\tau `$ has a branching ratio $`Br>3.51\times 10^3`$.

To see the contributions of the radiative decays precisely we present the radiative decay branching ratios with a cut of the photon energy, i.e., the branching ratios of the radiative decays $`D_sl\nu \gamma `$ with the photon energy $`E_\gamma k_{min}`$ as the follows: $`k_{min}=0.00001`$ GeV, $`k_{min}=0.0001`$ GeV, $`k_{min}=0.001`$ GeV, $`k_{min}=0.01`$ GeV and $`k_{min}=0.1`$ GeV respectively in Table (2). We also show the existing results of other methods in the same table.

Table (2): The Radiative Decay branching ratios with cuts of the photon momentum and the results of Ref

| $`k_{min}`$ | $`Br_e`$ | $`Br_\mu `$ | $`Br_\tau `$ |
| --- | --- | --- | --- |
| GeV | $`10^5`$ | $`10^4`$ | $`10^6`$ |
| 0.00001 | 2.552 | 4.901 | 6.336 |
| 0.0001 | 2.552 | 3.908 | 4.597 |
| 0.001 | 2.551 | 2.915 | 2.868 |
| 0.01 | 2.549 | 1.927 | 1.217 |
| 0.1 | 2.475 | 0.971 | 0.727 |
| Ref | 10 | 1 | $``$ |
| Ref | 17 | 1.7 | $``$ |
| Ref | 7.7 | 2.6 | 320 |

Considering the possible experimental resolutions of photon, our prediction of the radiative decay branching ratios $`Br(D_se\nu _e\gamma )`$ and $`Br(D_s\mu \nu _\mu \gamma )`$ are close to the values in Ref, but our prediction of $`Br(D_s\tau \nu _\tau \gamma )`$ is much smaller than the one in Ref. In our model, if we using a smaller cut $`k_{min}`$, then obtained a larger $`Br(D_s\mathrm{}\nu _{\mathrm{}}\gamma )`$, because the decay widths depend on $`Log(k_{min})\text{[8]}`$, the change of branching ratios will be not so much on the selection of $`k_{min}`$, we can see this in Table (2), and for another example, if $`k_{min}=1.0\times 10^{10}`$ GeV, then we obtain $`Br(D_se\nu _e\gamma )=3.59\times 10^5`$, $`Br(D_s\mu \nu _\mu \gamma )=1.38\times 10^3`$, $`Br(D_s\tau \nu _\tau \gamma )=2.11\times 10^5`$, but so small a $`k_{min}`$, it is very difficult in experiment. We can conclude that the best radiative decay channel is easy to search in experiment is $`D_s\mu \nu _\mu \gamma `$.

For the convenience to compare with experiments, we present the photon spectrum of the radiative decays in Fig.2 and Fig.3. In addition, we should note that the widths are quite sensitive to the decay constant $`f_{D_s}`$, and are sensitive to the values of the quark masses $`m_s`$ and $`m_c`$.


format: text

text:
# A lattice study of the nucleon excited states with domain wall fermions

## I Introduction

An important challenge in lattice calculations is to reproduce the hadron mass spectrum from first principles in quantum chromodynamics (QCD). The latest lattice QCD calculations of the light-hadron mass spectrum in the quenched approximation agree with experimental values within about 5% . However, this success is mainly restricted to ground states. Indeed, results are scarcely available for the excited-state mass spectrum.

Another essential shortcoming of these calculations which use Wilson or Kogut-Susskind fermions is the absence of chiral symmetry at finite lattice spacing, in accord with the Nielsen-Ninomiya no-go theorem . At non-zero lattice spacing Wilson fermions explicitly break the full chiral symmetry of the continuum down to the vector sub-group, so only flavor symmetry is preserved. On the other hand, Kogut-Susskind fermions have only a single exact $`U(1)`$ axial symmetry, and flavor symmetry is completely broken . Of course, it is expected that in the continuum limit, which is difficult to achieve in practice, both actions recover the full chiral symmetry.

Several years ago Kaplan constructed a new type of lattice fermion  known as domain wall fermions, which were further developed by Shamir  and also by Narayanan and Neuberger . Especially, the former reformulated it for lattice QCD simulations. The key feature of domain wall fermions is that they utilize an extra fifth dimension to circumvent the Nielsen-Ninomiya no-go theorem and maintain chiral symmetry at non-zero lattice spacing. In practical simulations the extra dimension is finite, so the chiral symmetry is not exact. The symmetry breaking is very soft, however, since it is highly suppressed with the number of sites in the extra dimension, $`L_s`$. In other words, $`L_s`$ gives us a way to control the violation of chiral symmetry. Domain wall fermions also possess exact flavor symmetry for any value of $`L_s`$.

Quenched lattice QCD calculations with domain wall fermions have shown that good chiral properties are obtained for moderate sizes of the fifth dimension, $`L_s1016`$, if the lattice spacing is small enough ($`a<0.1`$ fm) . Recent studies by the RIKEN-BNL-Columbia-KEK collaboration  and the CP-PACS collaboration  have quantified in detail the explicit chiral symmetry breaking effects due to finite $`L_s`$. For low energy QCD, the results can be simply summarized: there is a unique additive quark mass, $`m_{\mathrm{res}}`$, appearing in the low energy QCD lagrangian which arises from the finite size of the extra dimension . This additive quark mass has been measured quite accurately, and for $`L_s=16`$ and quenched $`\beta =6.0`$, $`m_{\mathrm{res}}`$ 3% of the strange quark mass. In this paper, we apply domain wall fermions to baryon excited states, especially the spin one-half and isospin one-half negative-parity nucleon, $`N^{}(1535)`$, as a further test of domain wall fermions in the baryon sector. For masses which are $`O(1)`$ GeV in the chiral limit, we do not expect $`m_{\mathrm{res}}`$ to have a significant effect.

We are interested in a long-standing puzzle in the excited state spectrum of the nucleon. The first question addressed in this paper is whether the mass difference between the nucleon $`N(939)`$ and the negative-parity nucleon $`N^{}(1535)`$ is well reproduced in lattice QCD. The spin one-half $`N^{}`$ state can be considered the parity partner of the nucleon. Of particular interest is the large mass splitting between $`N`$ and $`N^{}`$. From the viewpoint of parity partners, these two states would be degenerate if the relevant 2-flavor chiral symmetry were exact and preserved by the vacuum . Of course, there is no proof that this large mass splitting comes directly from the spontaneous breaking of chiral symmetry. However, at least spontaneous chiral symmetry breaking is responsible for the absence of such parity doubling since the explicit breaking is quite small in the case of two flavors. In this sense, regardless of a model or a theory, chiral symmetry and its spontaneous breaking are important for reproducing precisely the mass splitting between parity-partner hadrons. In this paper we show that domain wall fermions accurately reproduce the large observed mass splitting (some of our results have been reported earlier ).

Conventional lattice fermion schemes have had difficulty in this challenge. An early calculation  using Wilson fermions as well as a more recent one , both just managed to extract a mass splitting between the parity partners with large uncertainties.***In a recent preprint  which appeared after this work was complete it is clear that the chiral symmetry of their fermions was significantly improved by simulating closer to the continuum limit than previous studies which used Wilson fermions and by adding the clover term.

The first sophisticated calculation using an improved Wilson fermion action and relatively heavy quarks resolved a definite mass splitting between the parity partners which, however, was about a factor of two smaller than the observed splitting. Recent results using improved Wilson fermions have confirmed the large mass splitting over a wide range of quark masses that we have found with domain wall fermions. We note that the leading lattice spacing errors that are removed from these improved calculations break chiral symmetry. Although Kogut-Susskind fermions have a remnant $`U(1)`$ axial symmetry, they cannot be used practically for the spin one-half $`N^{}`$ mass calculation due to flavor mixing. The reason is that Kogut-Susskind fermions have only discrete flavor symmetries belonging to a subgroup of $`SU(4)`$  which contains three irreducible representations, 8, 8 and 16 for baryon operators. Two appropriate representations 8 and 16, to which $`N^{}`$(1535) belongs, also contain the negative-parity $`\mathrm{\Lambda }`$ states $`\mathrm{\Lambda }`$(1405) and $`\mathrm{\Lambda }`$(1520) and the spin three-half negative-parity nucleon state $`N^{}`$(1520) . Thus, the study of the spin one-half $`N^{}`$ spectrum with Kogut-Susskind fermions always faces inevitable contamination from lower mass states.

It is also interesting to note that a non-relativistic quark model with the so-called color magnetic interaction  and the MIT bag model , both of which explicitly break chiral symmetry, have some difficulty reproducing the large mass splitting between $`N(939)`$ and $`N^{}(1535)`$ without adoption of less realistic model parameters. The non-relativistic quark model is based on a harmonic oscillator description of the orbital motion of constituent quarks. As remarked in the appendix of the original paper by Isgur and Karl , the plausible value of its oscillator quantum should be roughly 250 MeV to reproduce the observed charge radius and magnetic moment of the nucleon. Since this model regards $`N^{}(1535)`$ as a state with one quantum excitation in orbital motion, it indicates that the corresponding $`N^{}`$ state lies at most a few hundred MeV above the ground state. Even worse, there is the serious problem of the wrong ordering between $`N^{}(1535)`$ and the positive-parity excited nucleon $`N^{}(1440)`$ because the corresponding $`N^{}`$ state should be assigned two oscillator quanta in this model . This wrong ordering problem does not seem to be easily alleviated . It is easy to see that the MIT bag model faces essentially the same problem, as the single quark states in the model alternate in parity with roughly even spacings .

The question arises how does this ordering of $`N^{}`$ and $`N^{}`$ appear in lattice QCD calculations? Until  this question could not be answered for lack of systematic results and investigations. A calculation of the mass of the positive-parity excited nucleon is, of course, much more difficult than the nucleon ground-state. Attempts have been made to evaluate the $`N^{}`$ mass from a two state fit to the nucleon correlation function . However, large statistics are required compared to a single exponential fit. Also it is difficult to control the systematic errors. Leinweber used the QCD sum rule (QCDSR) continuum model and quenched lattice data to estimate the threshold of contributions from excited states . However, the scaling of his result has been questioned by Allton and Capitani . In this paper we take an alternative approach, using the continuum-like behavior of the domain wall fermion operator  to obtain the mass of the positive-parity excited nucleon. Our results have been confirmed in  where a similar approach is employed in the context of improved Wilson fermions.

In conventional lattice QCD calculations an interpolating operator which is strongly associated with the non-relativistic limit is used to extract the mass of the nucleon ground state. A second unconventional operator, which does not have a non-relativistic limit, is discarded since it is expected to couple weakly to the ground state. Indeed, using this operator, no one has succeeded in evaluating the mass of the nucleon ground state in lattice QCD calculations with Wilson fermions  because of its small coupling to the ground state . The expectation of an approximately zero overlap on the ground state provides the possibility that the use of the second operator in lattice calculations directly yields the mass of the $`N^{}`$ state, at least in the relatively heavy (valence) quark mass region. This prospect, however, is built on the assumption that the lattice defined operators inherit the features of the continuum ones. In the case of the Wilson fermions, the Wilson term, which explicitly breaks chiral symmetry, induces mixing between the conventional and unconventional operators . Thus, the desired feature of the unconventional operator in the continuum is diminished in lattice calculations with Wilson fermions. On the other hand, this type of mixing between three quark operators is absent at one loop in perturbation theory using domain wall fermions with large $`L_s`$ ; thus we expect this mixing to be suppressed in domain wall fermion lattice calculations. Indeed, we find that the second operator has negligible overlap with the ground state and furthermore provides a reliable signal of the excited state as the asymptotic state in the heavy quark mass region.

The organization of our paper is as follows. In Sec. II, we briefly review the basic formulae and notation regarding domain wall fermions. In Sec. III we investigate the properties of the two-point correlation function for the nucleon and its parity partner. Sec. IV gives the details of our Monte Carlo simulations and the results for the parity partner of the nucleon ($`N^{}`$) and the first positive-parity excited nucleon ($`N^{}`$). Then, we compare our results to the experimental values. Finally, we present our conclusions in Sec. V.

## II Domain wall fermions

In this section we closely follow the development of domain wall fermions by Shamir . The domain wall fermion action is essentially regarded as a five-dimensional extension of the Wilson fermion action:

$$S_{\mathrm{DWF}}=\underset{x,x^{}}{}\underset{s,s^{}}{}\overline{\mathrm{\Psi }}(x,s)[\delta _{s,s^{}}D_{x,x^{}}^{}+\delta _{x,x^{}}D_{s.s^{}}^{}]\mathrm{\Psi }(x^{},s^{}),$$
(1)
where $`x,x^{}`$ are four-dimensional Euclidean space-time coordinates and $`s,s^{}`$ denote coordinates in the extra dimension labeled from 0 to $`L_s1`$ (to take advantage of existing high-speed computer code, our domain wall fermion Dirac operator is the Hermitian conjugate of the one in , hence our notation is the same as in ). Here, $`D_{x,x^{}}^{}`$ corresponds to the four-dimensional Wilson-Dirac operator with a mass term (domain wall height) $`M_5`$.

$$D_{x,x^{}}^{}=\frac{1}{2}\underset{\mu }{}\left\{(1\gamma _\mu )U_\mu (x)\delta _{x+\widehat{\mu },x^{}}+(1+\gamma _\mu )U_\mu ^{}(x^{})\delta _{x\widehat{\mu },x^{}}\right\}+(M_54)\delta _{x,x^{}},$$
(2)
where the Wilson term and mass term have opposite relative sign compared to the conventional one. $`D_{s,s^{}}^{}`$ is the five-dimensional analog of the four dimensional Wilson hopping term with $`\gamma _\mu `$ replaced by $`\gamma _5`$ and $`U_5(x,s)=1`$.

$`D_{s,s^{}}^{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{(1\gamma _5)\delta _{s+1,s^{}}+(1+\gamma _5)\delta _{s1,s^{}}2\delta _{s,s^{}}\right\}`$ (4)
$`{\displaystyle \frac{m_f}{2}}\left\{(1\gamma _5)\delta _{s,L_s1}\delta _{0,s^{}}+(1+\gamma _5)\delta _{s,0}\delta _{L_s1,s^{}}\right\}.`$
Note that the five dimensional fermions $`\mathrm{\Psi }(x,s)`$ are coupled only to four-dimensional gauge fields. The boundaries $`s=0`$ and $`s=L_s1`$ are anti-periodic and coupled with a weight $`m_f`$.

For an appropriate choice of $`M_5`$, this action has two chiral zero-modes of opposite handedness, one localized on each boundary of the fifth dimension. To simulate low energy QCD, four dimensional quarks are interpolated from the chiral modes,

$`q(x)`$ $`=`$ $`P_L\mathrm{\Psi }(x,0)+P_R\mathrm{\Psi }(x,L_s1),`$ (5)
$`\overline{q}(x)`$ $`=`$ $`\overline{\mathrm{\Psi }}(x,L_s1)P_L+\overline{\mathrm{\Psi }}(x,0)P_R,`$ (6)
where $`P_{R,L}=(1\pm \gamma _5)/2`$. With this definition, one can see that the $`m_f`$ terms in $`D_{s,s^{}}^{}`$ directly yield the usual four dimensional explicit chiral symmetry breaking term, $`m_f\overline{q}(x)q(x)`$, on the four dimensional layers $`s=0`$ and $`s=L_s1`$. The above definition of the quark fields is the simplest one but is not unique.

In the free theory, for the choice $`0<M_5<2`$, the domain wall fermion action corresponds to one flavor, and the four-dimensional light quark mass is given as $`m_q=m_fM_5(2M_5)`$ in the limit $`L_s\mathrm{}`$ and $`m_f0`$. This situation is approximately unchanged in lattice QCD simulations if $`M_5`$ is simply shifted, $`M_5M_5M_c`$ , where a simple estimate of $`M_c`$ is given in terms of the critical hopping parameter for four-dimensional Wilson fermions .

In the case of domain wall fermions, the $`SU(N_f)`$ axial transformation can be defined vectorially  as

$`[𝒬_A^a,\mathrm{\Psi }(x,s)]`$ $`=`$ $`+iϵ(s)\lambda ^a\mathrm{\Psi }(x,s),`$ (7)
$`[𝒬_A^a,\overline{\mathrm{\Psi }}(x,s)]`$ $`=`$ $`iϵ(s)\overline{\mathrm{\Psi }}(x,s)\lambda ^a,`$ (8)
where $`\mathrm{\Psi }(x,s)`$ is a five-dimensional zero mode. Opposite axial charges are assigned to fermions in the two half-spaces, so $`ϵ(s)=1`$ for $`0s<L_s/2`$ and $`ϵ(s)=1`$ for $`L_s/2sL_s1`$. The reason for this is clear: left and right handed modes are globally separated in the extra dimension.

With this axial transformation on the five dimensional fields, the four dimensional quark fields obey the familiar axial transformation:

$`[𝒬_A^a,q(x)]`$ $`=`$ $`+i\gamma _5\lambda ^aq(x),`$ (9)
$`[𝒬_A^a,\overline{q}(x)]`$ $`=`$ $`i\overline{q}(x)\gamma _5\lambda ^a.`$ (10)
In general the domain wall fermion action is not invariant under the transformations (7) and (8) in the limit $`m_f=0`$ because of a non-vanishing divergence of the axial current on the intermediate layers $`s=L_s/21`$ and $`s=L_s/2`$ which gives rise to an extra pseudo-scalar density in the axial Ward-Takahashi identity . However, such an anomalous term, which provides a chiral symmetry breaking effect due to mixing of the left- and right-handed modes, vanishes as $`L_s\mathrm{}`$ .

Theoretically, this residual breaking effect can be described by an additive quark mass $`m_{\mathrm{res}}`$ in the four-dimensional low energy effective Lagrangian for QCD , and recent simulations, in which $`m_{\mathrm{res}}`$ was determined in several ways, appear to confirm this. Simulations also show that for sufficiently small lattice spacing, $`m_{\mathrm{res}}0`$ as $`L_s\mathrm{}`$, and at the very least, $`m_{\mathrm{res}}`$ is small and accurately known for a wide range of $`L_s`$. Thus, in this study we will henceforth ignore these small effects and assume the chiral limit is $`m_f=0`$ instead of $`m_f=m_{\mathrm{res}}`$. Finally, we note that the non-perturbative origin of these effects and their relation to certain non-perturbative gauge field configurations is a very interesting and active area of research .

## III Baryon spectrum

### A Baryon two-point correlator

The mass $`m_B`$ of the low-lying baryon $`B`$ can be extracted from the two-point correlation function composed of the baryon interpolating operator $`𝒪_B`$, which has the appropriate quantum numbers specified by the desired state. Let us consider the vacuum expectation value of the time-ordered product of interpolating operators. The Euclidean time correlation function is projected out at zero spatial momentum through the sum over $`\stackrel{}{x}`$:

$$G_{𝒪_B}(t)=\underset{\stackrel{}{x}}{}0|T\{𝒪_B(\stackrel{}{x},t)\overline{𝒪}_B(0,0)\}|0,$$
(11)
which is dominated by the contribution of the lowest mass state for large Euclidean time, $`G_{𝒪_B}(t)\mathrm{exp}(M_B|t|)`$. For a spin one-half baryon, of course, the correlation function has non-trivial Dirac structure which may be expressed in the form

$$G_{𝒪_B}(t)\underset{\mathrm{large}t}{}(1+\gamma _4)\theta (t)A_Be^{M_Bt}+(1\gamma _4)\theta (t)A_Be^{+M_Bt},$$
(12)
where the positive definite $`A_B`$ is proportional to the square of the coupling strength between the interpolating operator $`𝒪_B`$ and the lowest mass state. The first and second terms are particle and anti-particle contributions, respectively. It is easy to get only the particle contribution by taking the trace with the projection operator $`P_+=(1+\gamma _4)/2`$, which we abbreviate as $`𝒪_B(t)\overline{𝒪}_B(0)=\mathrm{Tr}\{P_+G_{𝒪_B}(t)\}`$ hereafter.

### B Interpolating operators

Let us focus on the nucleon channel specified by the spin one-half iso-doublet non-strange baryons. In order to project out the desired channels, we have to construct interpolating operators from the quark fields with the appropriate quantum numbers ($`J=1/2`$ and $`I=1/2`$). However, there is considerable freedom in choosing the specific form of the composite operators. Indeed, there are two possible interpolating operators for the $`J^P=1/2^+`$ state, even if we restrict them to contain no derivatives and to belong to the $`(\frac{1}{2},0)(0,\frac{1}{2})`$ chiral multiplet under $`SU(2)_LSU(2)_R`$ :

$`B_1^+(x)`$ $`=`$ $`\epsilon _{abc}[u_a^T(x)C\gamma _5d_b(x)]u_c(x),`$ (13)
$`B_2^+(x)`$ $`=`$ $`\epsilon _{abc}[u_a^T(x)Cd_b(x)]\gamma _5u_c(x),`$ (14)
where $`abc`$ and $`ud`$ have usual meanings as color and flavor indices. $`C`$ is the charge conjugation matrix and the superscript $`T`$ denotes transpose. Here, Dirac indices have been suppressed. In Eq.(13) and (14), the superscript “$`+`$” refers to positive parity since these operators transform as $`𝒫B_{1,2}^+(\stackrel{}{x},t)𝒫^{}=+\gamma _4B_{1,2}^+(\stackrel{}{x},t)`$ under parity. To be precise, the linear combinations $`B_1^+\pm B_2^+`$ belong to distinct $`(\frac{1}{2},0)(0,\frac{1}{2})`$ chiral multiplets under $`SU(2)_LSU(2)_R`$ .

The operator $`B_1^+`$ alone is usually used in lattice QCD calculations to extract the nucleon ground state since $`B_2^+`$ couples only weakly to the ground state due to its vanishing in the non-relativistic limit . In fact, nobody has succeeded in extracting the nucleon mass spectrum from the $`B_2^+`$ operator. In our calculation, we also confirm that the ground state cannot be extracted from $`B_2^+(t)\overline{B}_2^+(0)`$; however, we have had some success with respect to the excited-state mass spectrum which we discuss in Sec IV-C.

Multiplying the left hand side of the previous positive-parity operators by $`\gamma _5`$, we obtain the interpolating operators with negative parity, $`J^P=1/2^{}`$ :

$`B_1^{}(x)`$ $`=`$ $`\gamma _5B_1^+(x)=\epsilon _{abc}[u_a^T(x)C\gamma _5d_b(x)]\gamma _5u_c(x),`$ (15)
$`B_2^{}(x)`$ $`=`$ $`\gamma _5B_2^+(x)=\epsilon _{abc}[u_a^T(x)Cd_b(x)]u_c(x),`$ (16)
since $`𝒫B_{1,2}^{}(\stackrel{}{x},t)𝒫^{}=\gamma _4B_{1,2}^{}(\stackrel{}{x},t)`$. The important point to notice is the relation between the correlation functions of opposite parities,

$$G_{B^+}(t)=\gamma _5G_B^{}(t)\gamma _5,$$
(17)
since $`B_{1,2}^{}=\gamma _5B_{1,2}^+`$. Eq.(17) means that the two-point correlation function of the spin one-half baryon can couple to both positive- and negative-parity states. Thus, the general form of the two point function is

$`G_{B^+}(t)`$ $`=`$ $`(1+\gamma _4)\theta (t)A_{B^+}e^{M_{B^+}t}+(1\gamma _4)\theta (t)A_{B^+}e^{+M_{B^+}t}`$ (18)
$``$ $`(1+\gamma _4)\theta (t)A_B^{}e^{+M_B^{}t}(1\gamma _4)\theta (t)A_B^{}e^{M_B^{}t}.`$ (19)
Note that the backward propagating contributions correspond to the anti-particles of the forward propagating states with opposite parity. The desired state is obtained by choosing the appropriate projection operator, $`1\pm \gamma _4`$, and direction of propagation.

Next, consider a lattice with finite extent $`T`$ in the time direction and (anti-)periodic boundary conditions. Eq. 19 is replaced by

$`G_{B^+}(t)`$ $`=`$ $`(1+\gamma _4)A_{B^+}e^{M_{B^+}t}\pm (1\gamma _4)A_{B^+}e^{M_{B^+}(Tt)}`$ (20)
$``$ $`(1+\gamma _4)A_B^{}e^{M_B^{}(Tt)}(1\gamma _4)A_B^{}e^{M_B^{}t},`$ (21)
where the (lower) upper sign stands for (anti-)periodic boundary condition on $`0<t<T`$. The anti-particle of the opposite-parity state can propagate through the time direction boundaries, so one faces unwanted contamination from the opposite-parity state in extracting the mass of the desired state. This contamination is not unavoidable if a double exponential fit is used. However, such fits require very high statistics. This is not a serious issue in the measurement of the nucleon ground state since the contamination from the negative-parity nucleon is expected to be negligible due to the large mass splitting $`M_B^{}M_{B^+}`$. Of course, this same splitting does affect the extraction of $`M_B^{}`$. The problem is resolved by choosing appropriate boundary conditions in the time direction to prevent the wrap-around effect, or by increasing the time extent $`T`$ sufficiently and placing the interpolating operators far from the boundary. In this study we take the former approach. It is common to employ Dirichlet boundary conditions where link valuables in the time direction at $`t=0`$ and $`t=T1`$ are set to zero when calculating quark propagators (for example, see ). However, we use a linear combination of two quark propagators with periodic and anti-periodic boundary conditions in the time direction to produce forward propagating states. Although our approach requires two times as many fermion matrix inversions to calculate one quark propagator, it does not suffer from (unknown) reflection effects at the time boundaries induced by Dirichlet boundary conditions.

### C Chiral symmetry and parity doubling

At the end of this section, let us briefly review how unbroken chiral symmetry imposes parity doubling in the hadron spectrum . We will generalize the following argument in Appendix A. For the sake of simplicity, we consider a particular transformation of the $`SU(2)`$ chiral symmetry, $`[𝒬_A,u]=+i\gamma _5u`$ and $`[𝒬_A,d]=i\gamma _5d`$. Then, one can easily find that the two-point correlators $`B_1^\pm `$ and $`B_2^\pm `$ transform as

$$[𝒬_A,B_{1,2}^\pm (x)\overline{B}_{1,2}^\pm (0)]=i\{\gamma _5,B_{1,2}^\pm (x)\overline{B}_{1,2}^\pm (0)\}$$
(22)
in the chiral limit. Now suppose that the vacuum possesses chiral symmetry: $`𝒬_A|0=0`$. According to Eq.(22), the two-point correlation functions anti-commute with $`\gamma _5`$,

$$\{\gamma _5,G_{B^\pm }(t)\}=0.$$
(23)
Immediately, with the help of Eq.(17), one finds

$$G_{B^+}(t)=G_B^{}(t),$$
(24)
which means that parity doubling arises in the nucleon channel due to chiral symmetry . Of course, in the real world chiral symmetry is spontaneously broken, $`𝒬_A|00`$, so that such parity doubling does not occur in the actual spectrum . In this sense, the spontaneous breaking of chiral symmetry is responsible for the absence of parity doubling. Thus, it seems important to properly handle chiral symmetry and its spontaneous breaking in order to calculate precisely the mass splitting between the nucleon and its parity partner. Finally, it is important to note that above argument ignores possible consequences to the ’t Hooft anomaly condition .

## IV Numerical results

### A Computational details

We generate quenched QCD configurations on a $`16^3\times 32`$ lattice with the standard single-plaquette Wilson action at $`\beta =6/g^2=6.0`$. The quark propagator is calculated using domain wall fermions with a fifth dimension of $`L_s`$=16 sites and a domain wall height $`M_5`$=1.8. Additional details of the simulation can be found in . Quenched $`\beta =6.0`$ corresponds to a lattice cut-off of $`a^11.9`$ GeV from $`aM_\rho `$=0.400(8) in the chiral limit and spatial size $`La1.7`$ fm .

We work in Coulomb gauge and calculate quark propagators using wall sources and local sinks. We expect that the use of wall sources provides better overlap with the desired states. To extract the state with desired parity in the spin one-half baryon spectrum, we construct forward (backward) propagating quarks by taking the appropriate linear combination of propagators with periodic and anti-periodic boundary conditions in the time direction, as mentioned before. This procedure eliminates the backward (forward) propagating opposite parity state which wraps around the time boundary.

In addition, we use two sources for quark propagators on each configuration ($`t_{\mathrm{src}}`$ = 5 and $`t_{\mathrm{src}}^{}`$ = 27) to increase statistics. After the appropriate parity projections, the correlation functions of the corresponding nucleon state can be folded together (averaged as a function of distance from the respective source). Here, time-slices are labeled from 0 to 31. To be clear, we give the form of the nucleon two-point function with arbitrary source location $`t_{\mathrm{src}}`$ after eliminating the unwanted contributions across time boundaries:

$$G_{B^+}(tt_{\mathrm{src}})=\{\begin{array}{cc}(1+\gamma _4)A_{B^+}e^{M_{B^+}(tt_{\mathrm{src}})}(1\gamma _4)A_B^{}e^{M_B^{}(tt_{\mathrm{src}})}\hfill & (T>t>t_{\mathrm{src}}),\hfill \\ (1\gamma _4)A_{B^+}e^{M_{B^+}(t_{\mathrm{src}}t)}(1+\gamma _4)A_B^{}e^{M_B^{}(t_{\mathrm{src}}t)}\hfill & (t_{\mathrm{src}}>t>0),\hfill \end{array}$$
(25)
which is constructed with the positive-parity interpolating operator (either $`B_1^+`$ or $`B_2^+`$). In the time range $`T>t>t_{\mathrm{src}}`$, positive- and negative-parity nucleon states are propagating only forward in time. This means that there are only particle contributions in this region. Thus, $`P_\pm =(1\pm \gamma _4)/2`$ project out positive- and negative-parity states, respectively. Alternatively, for $`t_{\mathrm{src}}>t>0`$ only the anti-particles contribute to the correlation function. In this case, $`P_\pm `$ act in the opposite way to the former case. In either case, all we have to do is calculate either $`G_{B^+}`$ or $`G_B^{}`$ to extract the masses for both positive- and negative-parity states since we already know $`\mathrm{Tr}\{P_\pm G_{B^+}(t)\}=\mathrm{Tr}\{P_{}G_B^{}(t)\}`$ with the help of Eq.17. In fact, we verified this relation on each configuration and then used $`G_{B_1^+}`$ and $`G_{B_2^{}}`$ to extract masses.

In our analysis, we use correlation functions in the ranges $`t>5`$ and $`t<27`$, for the sources $`t_{\mathrm{src}}=5`$ and $`t_{\mathrm{src}}^{}=27`$, respectively. After folding the propagators together as described above, we use a single offset from the source, $`t`$, defined in the range $`0<t<27`$

We use 405 independent gauge configurations for the lightest two quark masses, $`m_f=0.02`$ and $`0.03`$, 305 configurations for the intermediate ones, $`m_f=0.04`$ and $`0.05`$, and 105 configurations for the heavier ones, $`m_f=0.0750.125`$. These bare quark masses correspond to mass ratios $`M_\pi /M_\rho 0.590.90`$ as shown in Table I. In this calculation, $`SU(2)`$-isospin symmetry is enforced by equating $`m_f=m_f^{(u)}=m_f^{(d)}`$, so that the flavor index will not be explicitly displayed hereafter. All calculations were done on the 600 Gflops QCDSP machine at the RIKEN BNL Research Center.

### B Parity partner of nucleon: $`N^{}`$

We first calculate the effective masses $`M_{\mathrm{eff}}`$ to find appropriate time ranges for fitting. The effective mass is defined by

$$M_{\mathrm{eff}}(t)=\mathrm{ln}\{C(t)/C(t+1)\},$$
(26)
where $`C(t)`$ stands for $`B_{1,2}^\pm (t)\overline{B}_{1,2}^\pm (0)`$. We look for a plateau, or time independent region, in this quantity to extract the ground state mass. For example, Figures 1(a)-1(c) show effective masses for the nucleon ($`B_1^+`$) and its parity partner ($`B_{1,2}^{}`$) at $`m_f=0.03`$, 0.05, and 0.10. In Figure 1, the effective mass plot shows a clear plateau for the $`N`$, and one that is not as good for the heavier $`N^{}`$. Statistical uncertainties in Figure 1 are estimated by a single elimination jack-knife method. The effective masses for the $`N^{}`$ from both $`B_1^{}(t)\overline{B}_1^{}(0)`$ and $`B_2^{}(t)\overline{B}_2^{}(0)`$ agree well with each other for all three quark masses, except for $`m_f=0.03`$ and $`t<6`$ where there is a small difference outside of statistical errors.

In Figure 1 (a), the effective mass for the $`N^{}`$ becomes so noisy after $`t=9`$ that the value of the two-point correlator is consistent with zero within large errors, so we have left these points off of the plot. In addition, the effective mass shows a steep rise with increasing uncertainty as $`t`$ increases. This rise in the $`N^{}`$ effective mass weakens at relatively heavy quark masses where the signal becomes stable over 12 time-slices.

Next we present mass estimates of the $`N`$ and $`N^{}`$ obtained from covariant single exponential fits to the corresponding correlators. We fit each correlator from some minimum time-slice, $`t_{\mathrm{min}}`$, to an appropriate maximum time-slice ($`t_{\mathrm{max}}=20`$ for the $`N`$ and $`t_{\mathrm{max}}=1015`$ for the $`N^{}`$). $`t_{\mathrm{max}}`$ is roughly fixed with reference to the effective mass calculation. To keep fitting ranges as wide as possible, $`t_{\mathrm{min}}`$ is reduced from $`t_{\mathrm{max}}2`$ until $`\chi ^2/N_{\mathrm{DF}}>1.5`$ where $`N_{\mathrm{DF}}`$ denotes the degrees of freedom in the fit. Fitting details are given in Tables II-IV. All of our fits have confidence-level larger than 0.2 and estimates from the weighted average of the effective mass agree with the fitted masses within errors. A summary of our $`N`$ and $`N^{}`$ masses is given in Tables II-IV.

We can roughly estimate the systematic error coming from the choice of fitting range by varying $`t_{\mathrm{min}}`$. Let us determine the difference between our final fits and fits where the smallest time slice is not included. In the case of $`B_1^+`$ (the nucleon), resulting errors are much smaller than the statistical errors of the final fits. On the other hand, for both $`B_1^{}`$ and $`B_2^{}`$, the estimated systematic errors are comparable to the statistical errors, while for $`B_1^{}`$ at $`m_f=0.04`$ and 0.05 the systematic error is a factor of two larger than the statistical one. This uncertainty is directly related to the rise of the $`N^{}`$ effective mass, or poor plateau, in those cases where the correlation function was not statistically well resolved at larger times. This, in turn, has forced us to extract masses from the smaller $`t`$ region which may be contaminated by excited states.

In Figure 2 we show the low-lying nucleon spectrum as a function of the quark mass, $`m_f`$. The nucleon mass is extracted from $`B_1^+`$. We omit the point at $`m_f=0.02`$ for the operator $`B_2^{}`$ since a good plateau in the effective mass plot is absent. The $`N^{}`$ mass estimates from $`B_1^{}(t)\overline{B}_1^{}(0)`$ and $`B_2^{}(t)\overline{B}_2^{}(0)`$ agree with each other within errors in the whole quark mass range, as expected from their common quantum numbers . Note that this result disagrees with that obtained in : we find no discrepancy between masses extracted from $`B_1^{}(t)\overline{B}_1^{}(0)`$ and $`B_2^{}(t)\overline{B}_2^{}(0)`$. However, recent results in  and  are in good agreement with ours. We also obtain the same mass for the $`N^{}`$ from a mixed correlator $`B_1^{}(t)\overline{B}_2^{}(0)+B_2^{}(t)\overline{B}_1^{}(0)`$.

The most remarkable feature in Figure 2, which was first reported in , is that the $`N`$-$`N^{}`$ mass splitting is observed over the whole range of quark mass values and grows as the quark mass is decreased. To illustrate this point clearly, we compare two mass ratios in Figure 3, one from the baryon parity partners $`M_N^{}/M_N`$ and the other from pseudo-scalar and vector mesons $`M_\pi /M_\rho `$. Experimental points  are marked with stars, corresponding to non-strange (left) and strange (right) sectors. In the strange sector we use $`\mathrm{\Sigma }^+`$ and $`\mathrm{\Sigma }(1750)`$ as baryon parity partners and $`K`$ and $`K^{}`$ for the mesons. The baryon mass ratio clearly grows with decreasing meson mass ratio, toward the experimental values . We did not include the charm sector ($`M_D/M_D^{}0.93`$) since the parity partner of $`\mathrm{\Sigma }_c^{++}(2455)`$ is not measured experimentally. On the other hand, from our results we estimate the mass of this state to be roughly 2.7 GeV (we have used a simple linear in the quark mass ansatz to extract the mass from our degenerate quark data).

Finally, we evaluate the $`N`$ and $`N^{}`$ masses in the chiral limit. Taking a simple linear extrapolation in the four lightest quark masses for $`B_1^+`$ and $`B_1^{}`$, we find $`M_N`$=0.57(1) and $`M_N^{}`$=0.85(5) in lattice units. Setting the scale from the calculated $`\rho `$-meson mass , we obtain $`M_N1.1`$ GeV and $`M_N^{}1.6`$ GeV in the chiral limit. If we use the scale set by the calculated nucleon mass, we obtain $`M_N^{}1.5`$ GeV. Either way the $`N^{}`$ mass is in good agreement with the experimental values within about 5-10 %. The above errors do not include systematic uncertainties due to finite volume, non-zero lattice spacing, and quenching effects. Studies of such systematic errors will be addressed in future calculations.

Recently it was suggested  that the $`N^{}`$ propagator in quenched QCD may exhibit non-analytic chiral behavior associated with the anomalous contribution of the “$`\eta ^{}`$-$`N`$” intermediate state. Recall that in the quenched approximation the $`\eta ^{}`$ is also a pseudo-Goldstone boson. We stress that the quenched “$`\eta ^{}`$-$`N`$” intermediate state would give a negative metric contribution corresponding to a “quenched chiral loop” artifact . Needless to say, this spurious intermediate state, which is a source of unitarity violation, is not associated with the physical decay process .  This suggestion is inspired by recent articles . However, non-analytic effects which arise in the vicinity of the massless quark limit are hard to detect in our data since our calculation is not very close to this limit. In addition, such subtle effects could easily be mimicked by finite volume effects. The volume used in our study is still rather small. Furthermore, the statistical errors in our study are large enough, especially for the $`N^{}`$, that complicated fits beyond a naive linear one would not yield meaningful parameters or $`\chi ^2`$.

### C Unconventional nucleon operator

As mentioned earlier, the unconventional operator $`B_2^+`$ vanishes in the non-relativistic limit, so one may expect that it couples with negligible weight to the nucleon ground state near this limit. Indeed, no one has succeeded in extracting the ground state mass signal from $`B_2^+(t)\overline{B}_2^+(0)`$ in lattice QCD . Then, we expect that an approximately zero overlap on the ground state provides the possibility of direct access to the excited state using the $`B_2^+`$ operator. This prospect should hold as long as the lattice operator has the same symmetries as the continuum one.

In the massless quark limit, the combinations $`B_1^+\pm B_2^+`$ do not mix due to different chiral structure. In perturbation theory, each is multiplicatively renormalized with the same renormalization factor so that $`B_1^+`$ and $`B_2^+`$ also do not mix. However, conventional lattice fermions give rise to mixing through explicit chiral symmetry breaking . Thus, $`B_2^+`$ will couple to the ground state through unwanted mixing with $`B_1^+`$. On the other hand, the explicit breaking of chiral symmetry in domain wall fermions is highly suppressed, so there is hope that we may extract the positive-parity excited state of the nucleon from $`B_2^+`$ as the asymptotic state.

Let us first compare the effective mass plots of $`B_1^+(t)\overline{B}_1^+(0)`$ and $`B_2^+(t)\overline{B}_2^+(0)`$ correlators (see Figures 4a-4c). The correlator $`B_2^+(t)\overline{B}_2^+(0)`$ is considerably noisier so that only time-slices near the source are useful. Nevertheless, the effective mass from $`B_2^+(t)\overline{B}_2^+(0)`$ yields a plateau albeit with large statistical errors. The plateau becomes more satisfactory for heavier quark mass. These plateaus are obviously different from those extracted from the $`B_1^+`$ correlator. When we apply the single exponential fit to the two-point correlation function composed of $`B_2^+`$, we obtain a mass that is quite large compared to the mass extracted from $`B_1^+(t)\overline{B}_1^+(0)`$. In  similar results from the $`B_2^+`$ correlator have been found.

As mentioned earlier, an explanation of the above result  is that $`B_2^+`$ has negligible overlap with the nucleon ground state since it does not have a non-relativistic limit and thus provides a direct signal for the positive-parity excited state of the nucleon. Of course, this explanation is valid only in the heavy valence-quark mass limit.

We investigate further the possibility that $`0|B_2^+|N0`$ through the calculation of the mixed correlation function $`B_1^+(t)\overline{B}_2^+(0)+B_2^+(t)\overline{B}_1^+(0)`$. If it is true that the overlap with the nucleon $`0|B_2^+|N`$ becomes small with increasing valence-quark mass, mass estimates from the mixed correlator $`B_1^+(t)\overline{B}_2^+(0)+B_2^+(t)\overline{B}_1^+(0)`$ should be consistent with the nucleon for lighter quark mass and the positive-parity excited state for heavier quark mass. Such behavior is evident in Figure 5. We stress that the single particle fit to the mixed correlation function does not in general yield the mass of an asymptotic state, only in the limits discussed above, which is clear from the figure. We also note that for improved Wilson fermions the mixed correlation function yields the ground state mass for all quark masses, even heavy ones . Thus it appears in that case that there is still significant mixing of $`B_1^+`$ and $`B_2^+`$.

As a result, it is possible to identify $`B_2^+`$ with the positive-parity excited nucleon ($`N^{}`$) for heavy quarks (e.g. $`m_f>0.075`$ in our study) . Indeed, we see $`|0|B_2^+|N/0|B_2^+|N^{}|^210^3`$ from double exponential fits for $`B_2^+(t)\overline{B}_2^+(0)`$ at $`m_f=0.10`$ and 0.125. Of course, this feature weakens in the lighter quark mass region (from around $`m_f=0.05`$). However, even for $`m_f<0.075`$, mass estimates from $`B_2^+`$ are still considerably larger than the nucleon mass, so we cannot rule out the possibility that $`B_2^+`$ provides a signal for the positive-parity excited nucleon at even lighter quark mass values.

### D Diagonalization method for 2$`\times `$2 matrix correlator

Let us consider the $`N^{}`$ in more detail. The mass of the excited nucleon may be obtained from a two state fit to the $`B_1^+`$ correlator which, however, requires large statistics and that neither the ground state or the first excited state lies close to any other state. Attempts to extract the $`N^{}`$ mass using this method have failed to reproduce the observed mass . In our fits the Hessian matrix often becomes singular, so we cannot extract the excited state mass from double exponential fits. Instead, we take an alternative approach proposed in . First, we define the 2$`\times `$2 matrix correlator $`𝒞(t)`$ using the two distinct baryon operators

$$𝒞(t)=\left[\begin{array}{cc}c_{11}(t)& c_{12}(t)\\ c_{21}(t)& c_{22}(t)\end{array}\right],$$
(27)
where $`c_{ij}(t)=B_i^\pm (t)\overline{B}_j^\pm (0)`$. Next we write $`C(t)`$ in terms of a transfer matrix $`\lambda (t,t_0)`$,

$$𝒞(t)\psi =\lambda (t,t_0)𝒞(t_0)\psi ,$$
(28)
where $`t_0`$ is fixed and $`t>t_0`$. If only two states are propagating in a given system, the masses of the two states $`E_\alpha `$ ($`E_1>E_0`$) are given by the eigenvalues of the transfer matrix:

$$\lambda _\alpha (t,t_0)=e^{(tt_0)E_\alpha }(\alpha =0,1),$$
(29)
where $`E_\alpha `$ is independent of $`t_0`$. The smaller eigenvalue ($`\lambda _1`$) and larger eigenvalue ($`\lambda _0`$) refer to the masses of the excited state and the ground state respectively. In general, the system may have more than two states. Thus, we assume that two states become effectively dominant for an appropriately large time-slice $`t_0`$, which can be determined by checking the sensitivity of $`E_\alpha `$ with respect to variations of $`t_0`$.

We calculate the eigenvalues $`\lambda _\alpha (t,t_0)`$ of $`𝒞(t_0)^1𝒞(t)`$ for $`t_0=3`$, and then evaluate the effective masses of $`N`$ and $`N^{}`$ from the larger and smaller eigenvalues, respectively, with the statistical errors coming from the jack-knife method. Using these error estimates, we determine the masses in Table VI from weighted averages over the time-slice ranges listed there. Figure 6 shows that the results are quite consistent with the masses determined from single exponential fits to $`B_2^+(t)\overline{B}_2^+(0)`$ and $`B_1^+(t)\overline{B}_1^+(0)`$. We have checked that the variations of $`t_0`$ around $`t_0=3`$ do not significantly affect the effective masses. For $`m_f0.03`$, we unfortunately could not extract the mass of the first excited nucleon since a good plateau in the effective mass plot is absent. We note that for $`m_f=0.04`$ there is still a large splitting in the eigenvalues, indicating the overlap of $`B_2^+`$ with the nucleon remains small.

In contrast to the positive-parity state, the eigenvalues in the negative-parity case appear degenerate, as shown in Figure 7. This will be discussed in more detail in the next section.

As we have seen, the first excited state may die out so quickly that only a few time-slices are available for evaluating the mass, even if the excited state is well separated from the ground state. To circumvent this quick damping, simulations performed on an anisotropic lattice where the temporal lattice spacing is finer than the spatial one may be useful. In fact, a recent lattice study  has shown this to be an effective way to extract masses of nucleon excited states.

### E Comparison with experiment

In the physical spectrum, we have another negative-parity state, $`N^{}(1650)`$, which is just above the lowest state $`N^{}(1535)`$. Although these states are quite close to each other, they are easily distinguished due to a peculiar decay mode of $`N^{}(1535)`$. It is well known that the decay rate of $`N^{}(1535)N+\eta `$ is comparable with $`N^{}(1535)N+\pi `$ even though the $`N\pi `$ decay mode is kinematically favored over $`N\eta `$. For the case of $`N^{}(1650)`$ the $`N\pi `$ decay is dominant; the branching ratio of $`N\eta `$ is only a few percent. Thus, $`N^{}(1535)N+\pi `$ seems to be anomalously suppressed. However, the corresponding states cannot be easily distinguished without knowledge of decay patterns if they lie close to each other.

As mentioned, for the negative-parity case our data shows that both eigenvalues of the transfer matrix are the same within statistical errors (see Figure 7). Of course, this does not rule out the possibility that the splitting might become clear with more statistics and in the lighter quark mass region. That is because the central value of the difference between the two eigenvalues increases in the lighter quark mass region. Furthermore, the value of about 100 MeV at $`m_f=0.03`$ seems to be large enough to reproduce experimental splitting in the chiral limit.

Nevertheless, unlike the positive-parity case, we have no indication of the presence of two independent negative-parity states in our calculation where the mixed correlator yields a consistent mass with both $`B_1^{}`$ and $`B_2^{}`$ interpolating operators, as listed in Table VI. Thus, our combined results from all data for the negative-parity nucleon allows two other possibilities regarding $`N^{}(1650)`$. One is that $`N^{}(1535)`$ and $`N^{}(1650)`$ are completely degenerate or not independent in the quenched calculation. This situation resembles naive quark models and the MIT bag model where $`N^{}(1535)`$ and $`N^{}(1650)`$ are degenerate if we neglect the spin dependent interaction. In this case, the analysis through the diagonalization of the 2$`\times `$2 correlator is no longer helpful.

A Second possibility is that $`N^{}(1650)`$ is simply missing in our calculation, i.e. both $`B_1^{}`$ and $`B_2^{}`$ operators may not couple, or couple weakly, to the $`N^{}(1650)`$ state. This may be related to an argument in  where the authors show that the desired vanishing of the phenomenological $`\pi NN^{}`$ coupling for the $`N^{}(1535)`$ state results simply from the chiral transformation properties of particular interpolating operators $`B_1^{}`$ and $`B_2^{}`$ with the help of the soft pion limit. If their argument is relevant to the suppression of the $`N\pi `$ decay, it is possible that $`B_1^{}`$ and $`B_2^{}`$ couple strongly to $`N^{}(1535)`$ but not to $`N^{}(1650)`$. However, their argument is no longer applicable to a third interpolating operator

$$B_{3\mu }^{}=(u^TC\sigma _{\alpha \beta }u)\sigma _{\alpha \beta }\gamma _\mu d(u^TC\sigma _{\alpha \beta }d)\sigma _{\alpha \beta }\gamma _\mu u,$$
(30)
which belongs to the chiral multiplet $`(\frac{1}{2},1)(1,\frac{1}{2})`$ and has no derivative . The $`B_{3\mu }^{}`$ couples to both $`J=3/2`$ and $`J=1/2`$ states. In this case, $`N^{}(1650)`$ might be extracted from the $`B_{3\mu }^{}`$ operator since there is no reason for lack of coupling to the $`N^{}(1650)`$ state.

Finally, we would like to mention a remaining puzzle. We have done the first successful lattice calculation of both $`N^{}`$ and $`N^{}`$ spectra. As for the $`N^{}`$, we have reliable data only for relatively large values of $`m_f`$. Comparing the $`N^{}`$ mass with the $`N^{}`$ mass, we find that the ordering of $`N^{}`$ and $`N^{}`$ is inverted compared to experiment. Furthermore, the level spacing between $`N`$-$`N^{}`$ and $`N^{}`$-$`N^{}`$ is almost even . What we see here closely resembles the wrong ordering problem of the excited nucleon spectrum in naive quark models and the MIT bag model. However, if this result is true for heavy quarks, we can make an important prediction: the first excited state of the spin one-half $`\mathrm{\Sigma }_c`$ is a negative-parity state rather than a positive-parity state.

Unfortunately, our statistics did not allow the computation of the $`N^{}`$ mass in the light-quark region ($`m_f<0.03`$). In addition, our data show no evidence for the possibility of switching the level ordering between $`N^{}`$ and $`N^{}`$ towards the chiral limit. Indeed, we did not observe that the $`N^{}`$ mass decreased faster than the $`N^{}`$ mass with decreasing quark mass. However, finite volume effects for higher excited states may be more serious than for lower excited states. According to naive quark models or the MIT bag model, the $`N^{}`$ state is radially excited in contrast to the nucleon ground state and the $`N^{}`$ state. We note that the trend in the three heaviest quark mass points is at least consistent with the switching of the ordering in the chiral limit. Needless to say, the quenched approximation may also play a role in this puzzle. In summary, our results do not rule out the possibility of switching the ordering between $`N^{}`$ and $`N^{}`$ near the chiral limit, and to solve this remaining puzzle we need further systematic calculations toward the chiral limit.

## V Conclusion

We have studied the mass spectrum of the parity partner of the nucleon in quenched lattice QCD using domain wall fermions which preserve chiral symmetry to a high degree in lattice simulations. Most importantly we demonstrated that this method is capable of calculating the mass of the negative-parity $`N^{}`$ state in the spin one-half and isospin one-half baryon sector.

We made a systematic investigation of the $`N^{}`$ spectrum by using two distinct interpolating operators, $`B_1^{}`$ and $`B_2^{}`$. We found the $`N^{}`$ mass estimates extracted from them agree with each other. In practice the $`B_1^{}`$ correlator is more convenient than $`B_2^{}`$ in extracting the $`N^{}`$ mass because it is less noisy, especially in the light quark mass region.

We have found a definite mass splitting between $`N`$ and $`N^{}`$ states in the whole quark mass range we studied. Furthermore, this splitting grows with decreasing quark mass. The $`N^{}`$ mass and the $`N`$-$`N^{}`$ mass splitting in the chiral limit obtained by extrapolation are consistent with the experimental value within about 5-10%, depending on the mass we choose to set the lattice spacing, the nucleon or the $`\rho `$ meson. These results have been confirmed by subsequent lattice calculations using improved Wilson fermions . Needless to say this is very encouraging for further investigations of $`N^{}`$ physics using lattice QCD simulations.

In contrast to the negative parity operators, the positive parity operators, $`B_1^+`$ and $`B_2^+`$, yield distinct mass signals. From the $`B_1^+`$ operator we obtained a clean signal for the nucleon ground state, while from $`B_2^+`$ we always found a heavier mass. This is probably because the latter vanishes in the non-relativistic limit and thus has small overlap with the ground state. Indeed, we confirmed numerically that $`0|B_2^+|N0`$ and that $`B_2^+`$ yields the mass signal of the positive-parity excited state in the heavy quark region. We found that in the heavy quark mass region the mass of this excited state is about twice as high above the ground state mass as the $`N^{}`$ mass. This property of the $`B_2^+`$ operator was confirmed by comparing the results with the excited-state spectrum obtained from the diagonalization of a 2$`\times `$2 correlation matrix constructed from both interpolating operators and also by examining single particle fits to the mixed correlation function.

We did not resolve a long-standing puzzle regarding the excited-state spectrum of the nucleon, the inversion of the positive- and negative- parity first excited states with respect to experiment. However, we believe that the level switching between $`N^{}`$ and $`N^{}`$ must occur close to the chiral limit, though there is scant evidence in our results that this might happen. Taking our calculation at face value leads to the prediction that the first excited state of $`\mathrm{\Sigma }_c`$ has $`J^P=1/2^{}`$.

Needless to say, more work needs to be accomplished in order to achieve a fully systematic calculation of the excited nucleon spectrum. In the near future we plan to increase our volume and statistics to further explore the chiral limit, and to check that lattice spacing errors are small by running at several gauge couplings. We expect that the latter effect is small due to the chiral symmetry properties of domain wall fermions. A long term goal is to include dynamical fermions.

## Acknowledgment

The authors acknowledge useful discussions with all the members of the RIKEN-BNL-Columbia-KEK collaboration, but especially N. Christ, R. Mawhinney and A. Soni. All lattice simulations in this paper were done on the RIKEN BNL QCDSP supercomputer as part of the RIKEN-BNL-Columbia-KEK collaboration. We thank RIKEN, Brookhaven National Laboratory, and the U.S. Department of Energy for providing the facilities essential for the completion of this work. Finally, S.S. thanks all of the members of the RIKEN-BNL Research Center for their warm hospitality during his stay, where most of the present study has been carried out.

## Appendix A: General chiral transformation and parity partners

To extend our discussion of Sec.III-C to general chiral transformations, we prepare iso-doublet operators for $`B_1^+`$ and $`B_2^+`$,

$`_1^+(x)`$ $`=`$ $`\left(\begin{array}{c}\hfill \epsilon _{abc}[u_a^T(x)C\gamma _5d_b(x)]u_c(x)\\ \hfill \epsilon _{abc}[d_a^T(x)C\gamma _5u_b(x)]d_c(x)\end{array}\right),`$ (33)
$`_2^+(x)`$ $`=`$ $`\left(\begin{array}{c}\hfill \epsilon _{abc}[u_a^T(x)Cd_b(x)]\gamma _5u_c(x)\\ \hfill \epsilon _{abc}[d_a^T(x)Cu_b(x)]\gamma _5d_c(x)\end{array}\right).`$ (36)
The upper component corresponds to the proton and the lower component corresponds to the neutron. Also we can define iso-doublet operators for the negative parity nucleons as $`_{1,2}^{}=\gamma _5_{1,2}^+`$, respectively. Here we consider a general transformation of $`SU(2)_V`$ and $`SU(2)_A`$ symmetry,

$`[𝒬_V^a,q(x)]`$ $`=`$ $`i\tau _aq(x),`$ (37)
$`[𝒬_A^a,q(x)]`$ $`=`$ $`i\tau _a\gamma _5q(x),`$ (38)
where $`q=(u,d)^T`$ and $`\tau _a`$ is a Pauli matrix. Under these transformations, the two-point correlator composed of either $`_1^\pm `$ or $`_2^\pm `$ should transform as

$`[𝒬_V^a,_{1,2}^\pm (x)\overline{}_{1,2}^\pm (0)]`$ $`=`$ $`i[\tau _a,_{1,2}^\pm (x)\overline{}_{1,2}^\pm (0)],`$ (39)
$`[𝒬_A^a,_{1,2}^\pm (x)\overline{}_{1,2}^\pm (0)]`$ $`=`$ $`i\{\gamma _5\tau _a,_{1,2}^\pm (x)\overline{}_{1,2}^\pm (0)\}.`$ (40)
Eq.(40) along with the fact that $`𝒬_V^a|0=0`$ tells us that the nucleon two-point correlation function should be commute with $`\tau _a`$,

$$[\tau _a,G_^\pm (t)]=0,$$
(41)
which means that $`G_^+(t)`$ has only diagonal elements of iso-spin indices. In other words, the proton state and the neutron state are eigenstates of iso-spin, of course. If $`𝒬_A^a|0=0`$, it turns out that

$$\{\gamma _5\tau _a,G_^\pm (t)\}=\tau _a\{\gamma _5,G_^\pm (t)\}=0.$$
(42)
Thus, we can obtain the strict relation in terms of two-point correlators between the nucleon and its parity partner

$$G_^+(t)=G_{^{}}(t)$$
(43)
in the non-broken phase of chiral symmetry.

FIGURE CAPTIONS

The effective mass of the nucleon from the $`B_1^+(t)\overline{B}_1^+(0)`$ correlator ($`\times `$) and its parity partner from $`B_1^{}(t)\overline{B}_1^{}(0)`$ ($`\mathrm{}`$) and $`B_2^{}(t)\overline{B}_2^{}(0)`$ ($`\mathrm{}`$) correlators on an ensemble of 405 configurations at $`m_f=0.03`$. The solid lines and dashed lines represent each fitted mass and its statistical error.

The effective mass of the nucleon from the $`B_1^+(t)\overline{B}_1^+(0)`$ correlator ($`\times `$) and its parity partner from $`B_1^{}(t)\overline{B}_1^{}(0)`$ ($`\mathrm{}`$) and $`B_2^{}(t)\overline{B}_2^{}(0)`$ ($`\mathrm{}`$) correlators on an ensemble of 305 configurations at $`m_f=0.05`$. The solid lines and dashed lines represent each fitted mass and its statistical error.

The effective mass of the nucleon from the $`B_1^+(t)\overline{B}_1^+(0)`$ correlator ($`\times `$) and its parity partner from $`B_1^{}(t)\overline{B}_1^{}(0)`$ ($`\mathrm{}`$) and $`B_2^{}(t)\overline{B}_2^{}(0)`$ ($`\mathrm{}`$) correlators on an ensemble of 105 configurations at $`m_f=0.10`$. The solid lines and dashed lines represent each fitted mass and its statistical error.

$`N`$ ($`\times `$) and $`N^{}`$ ($`\mathrm{}`$ and $`\mathrm{}`$) masses versus the quark mass $`m_f`$ in lattice units ($`a^11.9`$GeV from $`aM_\rho `$=0.400(8) in the chiral limit). The corresponding experimental values for $`N`$ and $`N^{}`$ are marked with lower and upper stars. The $`N`$-$`N^{}`$ mass splitting is clearly observed. Symbols ($`\times `$, $`\mathrm{}`$ and $`\mathrm{}`$) are defined as in Figure 1.

Mass ratio of the negative-parity excited-state and positive-parity ground-state baryons versus mass ratio of the pseudoscalar meson and vector meson. All calculations are done for three degenerate valence-quarks. Ratios are calculated using fitted masses from the $`B_1^{}(t)\overline{B}_1^{}(0)`$ ($`\mathrm{}`$) and $`B_2^{}(t)\overline{B}_2^{}(0)`$ ($`\mathrm{}`$) correlators. Experimental points are denoted by stars.

The effective mass from the $`B_1^+(t)\overline{B}_1^+(0)`$ ($`\times `$) and $`B_2^+(t)\overline{B}_2^+(0)`$ ($``$) correlators on an ensemble of 405 configurations at $`m_f=0.03`$. The solid lines and dashed lines represent each fitted mass and its statistical error.

The effective mass from the $`B_1^+(t)\overline{B}_1^+(0)`$ ($`\times `$) and $`B_2^+(t)\overline{B}_2^+(0)`$ ($``$) correlators on an ensemble of 305 configurations at $`m_f=0.05`$. The solid lines and dashed lines represent each fitted mass and its statistical error.

The effective mass from the $`B_1^+(t)\overline{B}_1^+(0)`$ ($`\times `$) and $`B_2^+(t)\overline{B}_2^+(0)`$ ($``$) correlators on an ensemble of 105 configurations at $`m_f=0.10`$. The solid lines and dashed lines represent each fitted mass and its statistical error.

The fitted mass from $`B_1^+(t)\overline{B}_1^+(0)`$ ($`\times `$), $`B_2^+(t)\overline{B}_2^+(0)`$ ($``$) and mixed type $`B_1^+(t)\overline{B}_2^+(0)+B_2^+(t)\overline{B}_1^+(0)`$ ($`\mathrm{}`$) correlators. The corresponding experimental values for $`N`$ and $`N^{}`$ are marked with lower and upper stars. Note, the values extracted from the mixed correlation function do not represent the mass of an actual asymptotic state, except possibly in the light and heavy quark limits.

Comparison of the fitted mass from $`B_2^+(t)\overline{B}_2^+(0)`$ ($``$) and the estimated mass from the average effective mass of the smaller eigenvalue of the transfer matrix ($``$). The symbol $`\times `$ corresponds to the nucleon ground-state mass evaluated from the larger eigenvalue of the transfer matrix, which is quite consistent with the fitted mass from $`B_1^+(t)\overline{B}_1^+(0)`$.

Symbols $`\mathrm{}`$ and $`\mathrm{}`$ correspond to the the estimated mass from the average effective mass of the larger and smaller eigenvalues of the transfer matrix for the negative parity state. The eigenvalues are degenerate within errors. Symbols $`\times `$ and $``$ are defined as in Figure 6. The corresponding experimental values for $`N(939)`$, $`N^{}(1440)`$ and $`N^{}(1535)`$ are marked with lower, middle and upper stars. The ordering of the negative-parity nucleon ($`N^{}`$) and the positive-parity excited nucleon ($`N^{}`$) is inverted relative to experiment.

FIG.1a (Phys.Rev.D) Shoichi Sasaki et al.

FIG.1b (Phys.Rev.D) Shoichi Sasaki et al.

FIG.1c (Phys.Rev.D) Shoichi Sasaki et al.

FIG.2 (Phys.Rev.D) Shoichi Sasaki et al.

FIG.3 (Phys.Rev.D) Shoichi Sasaki et al.

FIG.4a (Phys.Rev.D) Shoichi Sasaki et al.

FIG.4b (Phys.Rev.D) Shoichi Sasaki et al.

FIG.4c (Phys.Rev.D) Shoichi Sasaki et al.

FIG.5 (Phys.Rev.D) Shoichi Sasaki et al.

FIG.6 (Phys.Rev.D) Shoichi Sasaki et al.

FIG.7 (Phys.Rev.D) Shoichi Sasaki et al.


text:
# Primary Sequences of Protein-Like Copolymers: Levy Flight Type Long Range Correlations

## Abstract

We consider the statistical properties of primary sequences of two-letter $`HP`$ copolymers ($`H`$ for hydrophobic and $`P`$ for polar) designed to have water soluble globular conformations with $`H`$ monomers shielded from water inside the shell of $`P`$ monomers. We show, both by computer simulations and by exact analytical calculation, that for large globules and flexible polymers such sequences exhibit long-range correlations which can be described by Levy-flight statistics.

PACS numbers: 61.41.+e, 36.20.Ey, 87.15.Cc

In refs , a new approach to the design of specific primary sequences for the $`HP`$-copolymers consisting of monomeric units of two types (hydrophobic $`H`$ and polar $`P`$) has been proposed by some of the present authors. Unlike some other methods of sequence design known in the literature (see review  and references therein), the approach in question does not aim to mimic folding into one particular conformation. The goal is to model simpler and more robust property of proteins, such as their ability to stay dissolved and shield their hydrophobic monomers from water. The essence of this approach is illustrated in Fig. 1. We start with an arbitrary computer-generated globular conformation of a homopolymer chain (formed due to the strong attraction of monomer units, Fig. 1a) and perform a ”coloring” procedure: monomer units in the core of the globule (having many neighbors) are set to be $`H`$-units while monomer units belonging to a globular surface (where the number of neighbors is smaller) are assigned to be of $`P`$-type, Fig. 1b. Then the obtained primary sequence is fixed, uniform attraction of monomer units is removed and newly generated $`HP`$-copolymer is ready for the further investigation (Fig. 1c). Thus obtained macromolecules are protein-like in the sense that they mimic segregation of globule into hydrophobic core and stabilizing hydrophilic envelope. The properties of protein-like copolymers were examined in ; see also  on possible ways of experimental realization.

In this Letter, we address correlations between $`H`$\- and $`P`$-units along the protein-like sequences. This may shed light on the conditions which must be met by the sequence to provide for the water solubility of globules - the issue of great potential relevance to our understanding of early evolution. We show, both by computer simulations and by exact analytical calculation, that correlations have a long-range character. More specifically, for the simple model of flexible polymer, they belong to the so-called Levy flight statistics.

To begin with, statistical properties of protein-like $`HP`$-sequences can be assessed computationally by the method similar to that used by Stanley et al  in their search for long-range correlations in DNA sequences. We choose the ”window” of length $`\mathrm{}`$, move it step by step along the generated $`HP`$-sequence, and at each step count the number of $`H`$ units inside the window. This number, which we write as $`_{i=j}^{j+\mathrm{}}u_i`$ is a random variable, depending on the position $`j`$ of the window along the sequence; here $`u_i`$ is the variable associated with every monomer $`i`$ such, that $`u_i=1`$ if monomer $`i`$ is $`H`$ and $`u_i=0`$ if it is $`P`$. This random variable has certain distribution. Its average is determined by the overall sequence composition (total numbers of $`H`$\- and $`P`$\- monomers), and its dispersion is easy to calculate:

$$D_{\mathrm{}}^2=\underset{i,j=k}{\overset{k+\mathrm{}}{}}\left(u_iu_ju_iu_j\right).$$
(1)
For a completely random $`HP`$-sequence, the value of $`D_{\mathrm{}}`$ scales as $`\mathrm{}^{1/2}`$ with the window width $`\mathrm{}`$. The dependence $`D_{\mathrm{}}\mathrm{}^\alpha `$ with $`\alpha >1/2`$ would then manifest the existence of long-range correlations.

The result of such calculation for averaging over 2000 independent protein-like $`HP`$-sequences of $`N=1024`$ monomer units with $`1:1`$ composition (obtained as in ref. ) is presented in Fig. 2, squares. For comparison, in the same figure 2 the data for two other types of sequences (averaged over 2000 independent species) are shown. One of them is purely random $`1:1`$ sequence, it demonstrates $`D_{\mathrm{}}\mathrm{}^{1/2}`$ scaling. Comparing this curve with Monte Carlo results we see immediately that protein-like sequence is not random and some correlations do exist in it. Thus, it is interesting to compare the squares in Fig. 2 and the dashed curve showing data for the sequence which was called ”random-block” in refs : the lengths of $`H`$ and $`P`$-blocks in a sequence are determined by the Poisson distributions adjusted to achieve the same $`1:1`$ composition and the same ”degree of blockiness” (average block length) as for protein-like $`HP`$-copolymer. This sequence exhibits a somewhat more rapid variation of $`D_{\mathrm{}}`$ at small $`\mathrm{}`$, but ultimately the law $`D_{\mathrm{}}\mathrm{}^{1/2}`$ is obeyed for large values of $`\mathrm{}`$. Nevertheless, this random-block model is also seen to be unsatisfactory for the statistical behavior of protein-like sequence throughout the interval of $`\mathrm{}`$ examined, $`2<\mathrm{}<500`$. Although the data do not fit accurately to any power law $`D_{\mathrm{}}\mathrm{}^\alpha `$, the slope of the observed $`D_{\mathrm{}}`$ dependence corresponds to $`\alpha `$ significantly larger than $`1/2`$, up to about $`0.85`$, thus indicating pronounced long-range correlations in protein-like sequence. In what follows we present analytical theory which produces curve $`D_{\mathrm{}}`$ in Fig. 2 in complete agreement with observations.

First of all, let us turn to the origin of long-range correlations in the primary sequences for $`HP`$-copolymers generated via the procedure illustrated in Fig. 1. Conceptually, this problem is fairly easy to address: since sequence in this scheme is uniquely determined by the parent conformation, the statistics of sequences reflects nothing but the statistics of parent conformations, which, in turn, is well understood. Indeed, the coloring procedure (Fig. 1b) operates in dense globular conformation. Since we consider very compact parent conformations, the statistics of polymer chain conformations inside the globule is ideal (Gaussian) according to the well-known Flory theorem . Therefore, all the statistical properties of parental conformations, including the correlations in the primary sequences produced by coloring procedure, can be derived via the solution of diffusion equation for random walks with appropriate boundary conditions.

To understand the fractal aspect of the sequences, it is convenient to concentrate on their uninterrupted homo-colored sections and on the points of connection between them. Our coloring procedure (Fig. 1) introduces a separation sphere of radius $`R^{}<R`$, such that all the units which in the parent conformation are confined inside this sphere are of $`H`$-type and

the units belonging to the shell layer $`R^{}<r<R`$ are of $`P`$-type. Therefore, homo-colored section is produced in our model by the chain section of the parent conformation placed entirely in either internal or shell regions of the globule. The probability to have an uninterrupted succession of some $`k`$ of $`H`$-monomer units in the sequence is equal to the probability that Gaussian (due to Flory theorem) polymer has a loop of $`k`$ monomer units entirely confined in $`H`$-region with ends on the separation surface. Similarly, the probability to have an uninterrupted succession of some $`k`$ of $`P`$-monomer units in the sequence is equal to the probability that ideal parental conformation has a loop of $`k`$ monomer units confined within the shell $`P`$-region, again with ends on the separation surface.

To address probability distributions $`P_H(k)`$ and $`P_P(k)`$, we begin with simple physical arguments yielding

$$P_{H,P}(k)\{\begin{array}{cc}k^{3/2},\hfill & 1<k<\left(\frac{d_{H,P}}{a}\right)^2\text{;}\hfill \\ \left(\frac{a}{d_{H,P}}\right)^3e^{\lambda _{H,P}k},\hfill & k>\left(\frac{d_{H,P}}{a}\right)^2\text{.}\hfill \end{array}$$
(2)
The upper asymptotic form is valid for short polymer loops, when neither curvature of the separating boundary nor overall globule shape play any role. In this regime, $`P(k)`$ is simply the probability for a random walk to start at the planar wall and to return to it for the first time after $`k`$ steps. This is classical probabilistic ”first return” problem, for which the $`k^{3/2}`$ answer is well known . This scaling is valid for loop sizes $`a\sqrt{k}`$ much larger than unity but smaller than the relevant characteristic length scale, $`d_H=R^{}`$ for the $`H`$-loops inside the inner sphere, or $`d_P=RR^{}`$ for the $`P`$-loops in the spherical shell. The second asymptotic form in equation (2) indicates that for long polymer loops the function $`P(k)`$ decays exponentially. It is easier to explain this in terms of polymer statistics: to confine a polymer chain of $`k`$ monomer units in a cavity costs some entropy $`\mathrm{\Delta }S`$, at $`a\sqrt{k}d`$ this entropy goes linearly with $`k`$, making the probability, $`\mathrm{exp}(\mathrm{\Delta }S)`$, exponential in $`k`$.

Let us now look closer at the cross-over values of $`k`$. In order to achieve the $`1:1`$ composition, $`R^{}`$ must be chosen such that volumes of internal $`H`$\- and shell $`P`$-regions are the same, which means $`R^{}=2^{1/3}R0.8R`$. The volume fraction of polymer units in a globule, $`\varphi `$, is controlled by the energy of interactions of monomer units used to prepare parental conformation. It is clear that $`R0.6aN^{1/3}\varphi ^1`$. Therefore, $`H`$-loops remain in the power law long range correlation regime up to the length $`k<0.24N^{2/3}\varphi ^2`$, while for $`P`$-loops this cross-over occurs somewhat earlier: $`k<0.015N^{2/3}\varphi ^2`$. Thus, we predict that there should be over a decade of length scales in which $`H`$-loops are still long range correlated, while only short range correlations remain in the $`P`$-loops.

The result (2) is sufficient to explain qualitatively correlations in protein-like sequences, including the data shown in Fig. 2. Indeed, according to our discussion, protein-like sequence can be thought of as an alternating succession of $`H`$\- and $`P`$-stretches, with lengths of stretches taken independently from the corresponding distributions $`P_H(k)`$ and $`P_P(k)`$. This mathematical scheme is called a Levy flight . We conclude that the long-range correlations in the primary sequences of protein-like copolymers are described by Levy flight statistics. Furthermore, for the $`k^{3/2}`$ behavior of $`P(k)`$, the averaged block length diverges, and, therefore, the value of $`D_{\mathrm{}}`$ in the power law regime is controlled by the longest block, yielding $`D_{\mathrm{}}\mathrm{}`$, or $`\alpha =1`$. This is true as long as both $`H`$\- and $`P`$\- loops remain in fractal regime. On the other hand, when all loops cross-over to exponential distribution, $`D_{\mathrm{}}`$ crosses over to $`\alpha =1/2`$:

$$D_{\mathrm{}}\{\begin{array}{cc}\mathrm{},\hfill & \text{for }1<\mathrm{}<0.015N^{2/3}\varphi ^2\hfill \\ \mathrm{}^{1/2},\hfill & \text{for }\mathrm{}>0.24N^{2/3}\varphi ^2\hfill \end{array}$$
(3)
The cross-over region for $`D_{\mathrm{}}`$ is very broad, it corresponds to the situation in which $`P`$ loops are already ”large, while $`H`$-loops are still ”small.” Both $`\alpha =1`$ and $`\alpha =1/2`$ limits and wide cross-over agree qualitatively well with computational data, Fig. 2. This motivates more careful theory, in which instead of scaling estimates (2) and (3), the expressions in terms of infinite series, suitable for numerical calculation, is obtained.

To develop full analytical theory, it is convenient to use the random walk terminology to describe parent conformation. In this language, for instance, $`P_H(k)`$ is the probability that the random walker enters a sphere of the radius $`R^{}`$ and then arrives back to the boundary for the first time after ”time” $`k`$. Recall that the statistical weight of all random walk trajectories starting at the point $`\stackrel{}{r}_0`$ and arriving after $`k`$ steps at the point $`\stackrel{}{r}`$, $`G(\stackrel{}{r},k|\stackrel{}{r}_0)`$, obeys the diffusion equation

$$\frac{G(\stackrel{}{r},k|\stackrel{}{r}_0)}{k}=\frac{a^2}{6}\mathrm{\Delta }G(\stackrel{}{r},k|\stackrel{}{r}_0)+\delta (k)\delta (\stackrel{}{r}\stackrel{}{r}_0)$$
(4)
where $`a^2`$ is the mean square length of one step (the squared size of one monomer unit along the chain). To introduce the condition of first return we have to say that the walker never touches the boundary, which is achieved by imposing the boundary condition

$$G(\stackrel{}{r},k|\stackrel{}{r}_0)|_{\left|\stackrel{}{r}\right|=R^{}}=0.$$
(5)
The probability distribution of the ”first return times” in terms of $`G`$ is then given as the time-dependent flux of diffusing particles through the absorbing wall:

$$P_H(k)=|d\sigma \frac{a^2}{6}\frac{G}{r}|_{r=R^{}}|,$$
(6)
where $`/r`$ means the component of gradient normal to the surface, integration is performed over the closed separating surface, and the absolute value is written to avoid thinking about the direction of the flux. The normalization condition $`P_H(k)𝑑k=1`$ is guaranteed by the fact that all diffusing particles eventually leave through the surface. As regards $`\stackrel{}{r}_0`$, it should be taken within a distance of order $`a`$ from the separating $`R^{}`$ surface. The problem thus formulated, including equations (4-6), is easy to solve: we write $`G`$ in terms of bilinear expansion $`G=_ne^{k\lambda _n}\psi _n(\stackrel{}{r})\psi _n(\stackrel{}{r}_0)`$ over the eigenfunctions $`\psi _n`$ satisfying $`(a^2/6)\mathrm{\Delta }\psi _n=\lambda _n\psi _n`$ with the boundary condition (5). Upon spherical integration in (6), all angular dependent harmonics vanish, and we arrive at

$`P_H(k)={\displaystyle \frac{\pi a^2}{3R^{}r_0}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}n(1)^{n+1}\mathrm{sin}\left(n\pi {\displaystyle \frac{r_0}{R^{}}}\right)`$ (7)
$`\times \mathrm{exp}\left[{\displaystyle \frac{a^2}{6}}\left({\displaystyle \frac{n\pi }{R^{}}}\right)^2k\right]`$ (8)
The distribution $`P_P(k)`$ can be derived similarly, except that now we have to take care of the boundary condition at the outer surface of the globule. To this end, we argue that this condition must be taken in the form

$$_\stackrel{}{r}G(\stackrel{}{r},k|\stackrel{}{r}_0)|_{r=R}=0$$
(9)
Indeed, formally this condition ensures the constant density of monomer units throughout the globule for large values of $`k`$, as well as breaking of correlations as soon as polymer chain is ”reflected” by a globular boundary. Physically, this boundary condition reflects the fact that there is always a ”sticky layer” (or depletion layer) formed self-consistently along the internal surface of the globule due to the effective attraction of monomer units to the outer region where polymer density is depleted and excluded volume effect is reduced. As long as we are not interested in the structure of surface layer of the globule, we can just replace this layer by the effective boundary condition (9). After calculations for $`P_P(k)`$ we obtain

$`P_P(k)={\displaystyle \frac{a^2R^{}}{3(RR^{})^2r_0}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\zeta _n^2\mathrm{sin}\left(\zeta _n\frac{r_0R^{}}{RR^{}}\right)}{\zeta _n\mathrm{sin}\zeta _n\mathrm{cos}\zeta _n}}`$ (10)
$`\times \mathrm{exp}\left[{\displaystyle \frac{a^2}{6}}\left({\displaystyle \frac{\zeta _n}{RR^{}}}\right)^2k\right],`$ (11)
where $`\zeta _n`$ satisfies $`\zeta _n=\left(1R^{}/R\right)\mathrm{tan}\zeta _n`$.

Finally, to compute the dispersion $`D_{\mathrm{}}`$, we note that $`u_iu_j`$ in eq. (1) is the probability that both units $`i`$ and $`j`$ are of $`H`$ type, which happens if both are located inside the $`R^{}`$ region in the parental conformation. Thus

$`u_iu_j`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle _V^{}}d^3r{\displaystyle _V^{}}d^3r_0G(\stackrel{}{r},\left|ij\right||\stackrel{}{r}_0)=`$ (12)
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}e^{\lambda _n\left|ij\right|}\left({\displaystyle _V^{}}d^3r\psi _n(\stackrel{}{r})\right)^2,`$ (13)
where $`G(\stackrel{}{r},k|\stackrel{}{r}_0)`$ is the Green function satisfying eq. (4) with the boundary condition (9), $`\psi _n`$ and $`\lambda _n`$ are the corresponding eigenfunctions and eigenvalues. Plugging this into the eq. (1), one arrives at the cumbersome looking expression for $`D_{\mathrm{}}`$ which is easy to implement numerically; the result is plotted in Fig. 2 and shows virtually perfect fit to the Monte Carlo data . In fact, this fit may even be somewhat fortuitous; indeed, along with the diffusion equation (4), which is an approximation for the underlying bond fluctuation model, we can also switch from summation to integration in eq. (1), yielding

$$D_{\mathrm{}}^2=6\mathrm{}^2\underset{n=1}{\overset{\mathrm{}}{}}g\left(\frac{\xi _n^2a^2}{6R^2}\mathrm{}\right)\alpha _n,$$
(14)
where $`\xi _n`$ satisfies the equation $`\xi _n=\mathrm{tan}\xi _n`$, $`\alpha _n=\frac{\xi _n^2+1}{\xi _n^6}\left[\frac{\xi _nR^{}}{R}\mathrm{cos}\frac{\xi _nR^{}}{R}\mathrm{sin}\frac{\xi _nR^{}}{R}\right]^2`$, and the function $`g`$ is defined as $`g(x)=2(x1+\mathrm{exp}(x))/x^2`$. (Note that the sum in eq. (14) starts from $`n=1`$ and does not include the ground state, for which $`\psi _0`$ is a constant). It is easy to check that equation (14) does indeed have asymptotic behavior in accord with (3), including a broad cross-over; similarly, eqs. (8,11) agree with (2). Numerically, eq. (14) fits pretty well to the data (Fig. 2), but, we repeat, best fit is achieved by the cumbersome discrete formula.

In conclusion, we have shown that protein-like copolymers generated according to the coloring procedure proposed earlier  exhibit long-range correlations in the primary sequences. For the flexible polymers and large enough globules, these correlations belong to the Levy-flight statistics. This result, first observed in computer experiments, is confirmed and explained by analytical calculation. Analytical theory suggests that Levy flight statistics, albeit with a broader cross-over region, is expected even if parental conformation is not maximally compact, but rather a globule somewhat closer to the $`\theta `$-point. It becomes clear from our model that segregation of globule into a hydrophobic core and a hydrophilic peel, which is the necessary condition for water solubility, does impose severe restrictions on the sequence, and, therefore, must be manifested in certain correlations. Qualitatively, this agrees with earlier study of protein sequences . However, precise form of correlations may be affected by both globule size and chain flexibility, including the aspect of secondary structure. Moreover, another effect exists which favors correlations of the opposite sign, and which dominates for small globules and/or rigid polymers . Identification of long range correlations in protein sequences becomes, therefore, an interesting task promising to shed light on the evolutionary criteria involved in the selection of proteins and the role of water solubility among them.

The work was supported by NATO (PST/CLG. 974956), INTAS (INTAS-OPEN-97-0678) and RFBR 98-03-33337a. The authors thank Dr. A.Irbäck for the useful correspondence.


text:
# Viscoelastic response of sonic band-gap materials

## Appendix A

Applying the proper boundary conditions at the interface between the surrounding fluid and the sphere, Brill  requiring the continuity of the radial component of the displacement field and the surface traction, along with the requirement that there is no tangential component of the surface traction at the interface, we can determine the scattering matrix which connects the incident with the scattered field. Here for reasons of completeness, along the lines of the formalism established in Ref. ipsarob1, , the nonzero elements of the T matrix for a solid sphere in a fluid host are

$$T_{lm;l^{}m^{}}^{LL}=\frac{W_l^{LL}}{D_l}\delta _{ll^{}}\delta _{mm^{}},l,l^{}0,$$
(6)
with $`z_l=Sq_l`$ referring to the fluid and $`x_\nu =Sq_{s\nu }`$, with $`\nu =l,t`$ to the sphere. The superscripts $`LL`$ refer to the case of scattering in a fluid host, since incident and scattered waves are $`L`$-type compressional waves. The $`3\times 3`$ determinants $`D_l`$ and $`W_l^{LL}`$ are given by

$`D_l=\left|\begin{array}{ccc}d_{22}\hfill & d_{23}\hfill & d_{24}\hfill \\ d_{32}\hfill & d_{33}\hfill & d_{34}\hfill \\ d_{42}\hfill & d_{43}\hfill & d_{44}\hfill \end{array}\right|,W_l^{LL}=\left|\begin{array}{ccc}d_2^L\hfill & d_{23}\hfill & d_{24}\hfill \\ d_3^L\hfill & d_{33}\hfill & d_{34}\hfill \\ d_4^L\hfill & d_{43}\hfill & d_{44}\hfill \end{array}\right|,`$ (13)
where

$$\begin{array}{c}d_{22}=z_lh_{l}^{+}{}_{}{}^{}(z_l),d_{23}=l(l+1)j_l(x_t),d_{24}=x_lj_l^{}(x_l),\hfill \\ d_{32}=0,d_{33}=\left[l(l+1)1x_t^2/2\right]j_l(x_t)x_tj_l^{}(x_t),\hfill \\ d_{34}=x_lj_l^{}(x_l)j_l(x_l),d_{42}=x_t^2\rho h_l^+(z_l)/(2\rho _s),\hfill \\ d_{43}=l(l+1)\left[x_tj_l^{}(x_t)j_l(x_t)\right],\hfill \\ d_{44}=\left[l(l+1)x_t^2/2\right]j_l(x_l)2x_lj_l^{}(x_l),\hfill \\ d_2^L=z_lj_l^{}(z_l),d_3^L=0,d_3^L=x_t^2\rho j_l(z_l)/(2\rho _s).\hfill \end{array}$$
(14)
$`j_l^{}`$ and $`h_{l}^{+}{}_{}{}^{}`$ denote the first derivatives of the spherical Bessel and Hankel functions, respectively. The $`𝐓`$ matrix, because of spherical symmetry is diagonal in $`l`$ and independent of $`m`$.

The exact form of the above $`𝐓`$ matrix allows us to compute the normalized total scattering cross section of an elastic sphere (scattering cross section over $`\pi S^2`$) in a fluid host,

$$\sigma _{total}=\frac{4}{z_l^2}\underset{l=0}{}(2l+1)\left|T_{lm;lm}^{LL}\right|^2.$$
(15)
In addition the normalized absorption cross section is given by

$$\sigma _{abs}=\frac{4}{z_l^2}\underset{l=0}{}(2l+1)\left\{\left|T_{lm;lm}^{LL}\right|^2+\mathrm{Re}[T_{lm;lm}^{LL}]\right\}.$$
(16)