text:
# Study of Time-Dependent 𝑪𝑷-Violating Asymmetries in 𝒃→𝒔𝒒̄𝒒 Decays

on leave from \] Fermi National Accelerator Laboratory, Batavia, Illinois 60510 on leave from \]University of Pittsburgh, Pittsburgh PA 15260 on leave from \]Nova Gorica Polytechnic, Nova Gorica The Belle Collaboration

## Abstract

We present a measurement of $`CP`$-violation parameters in the $`bs\overline{q}q`$ penguin transitions ($`q=s,u,d`$) based on a 78 fb<sup>-1</sup> data sample collected at the $`\mathrm{{\rm Y}}(4S)`$ resonance with the Belle detector at the KEKB energy-asymmetric $`e^+e^{}`$ collider. One neutral $`B`$ meson is reconstructed in the $`\varphi K_S^0`$, $`K^+K^{}K_S^0`$, or $`\eta ^{}K_S^0`$ decay channel, and the flavor of the accompanying $`B`$ meson is identified from its decay products. $`CP`$ violation parameters for each of the three modes are obtained from the asymmetries in the distributions of the proper-time intervals between the two $`B`$ decays.

preprint: KEK Preprint 2002-128preprint: Belle Preprint 2002-40

In the standard model (SM), $`CP`$ violation arises in weak interactions from an irreducible complex phase in the Kobayashi-Maskawa (KM) quark-mixing matrix bib:KM . In particular, the SM predicts $`CP`$-violating asymmetries in the time-dependent rates for $`B^0`$ and $`\overline{B}^0`$ decays to a common $`CP`$ eigenstate $`f_{CP}`$ bib:sanda . Recent measurements of the $`CP`$-violating parameter $`\mathrm{sin}2\varphi _1`$ by the Belle bib:CP1\_Belle ; bib:Belle\_sin2phi1\_78fb-1  and BaBar bib:CP1\_BaBar  collaborations established $`CP`$ violation in neutral $`B`$ meson decays mediated by the $`bc\overline{c}s`$ tree transition bib:CC  at a level consistent with KM expectations.

Despite this success, many tests remain before one can conclude that the KM model provides a complete description. For example, the charmless decays $`B^0\varphi K_S^0`$, $`B^0K^+K^{}K_S^0`$, and $`B^0\eta ^{}K_S^0`$, which are mediated by the $`bs\overline{s}s`$ transition ($`B^0\eta ^{}K_S^0`$ also receives $`bs\overline{d}d`$ and $`bs\overline{u}u`$ penguin contributions) are potentially sensitive to new $`CP`$-violating phases from physics beyond the SM bib:lucy . SM contributions from the $`bu\overline{u}s`$ tree diagram are expected to be highly suppressed bib:tree-penguin ; bib:Garmash . Thus, the SM predicts that $`CP`$ violation measurements in these charmless modes should yield $`\mathrm{sin}2\varphi _1`$ to a good approximation. Consequently, a significant deviation in the time-dependent $`CP`$ asymmetry in these modes from what is observed in $`bc\overline{c}s`$ decays would be evidence of a $`CP`$-violating phase not expected in the KM model.

In the decay chain $`\mathrm{{\rm Y}}(4S)B^0\overline{B}{}_{}{}^{0}f_{CP}f_{\mathrm{tag}}`$, where one of the $`B`$ mesons decays at time $`t_{CP}`$ to a final state $`f_{CP}`$ and the other decays at time $`t_{\mathrm{tag}}`$ to a final state $`f_{\mathrm{tag}}`$ that distinguishes between $`B^0`$ and $`\overline{B}^0`$, the decay rate has a time dependence given by bib:CPVrev

$$𝒫(\mathrm{\Delta }t)=\frac{e^{|\mathrm{\Delta }t|/\tau _{B^0}}}{4\tau _{B^0}}\left\{1+q\left[𝒮\mathrm{sin}(\mathrm{\Delta }m_d\mathrm{\Delta }t)+𝒜\mathrm{cos}(\mathrm{\Delta }m_d\mathrm{\Delta }t)\right]\right\},$$
(1)
where $`\tau _{B^0}`$ is the $`B^0`$ lifetime, $`\mathrm{\Delta }m_d`$ is the mass difference between the two $`B^0`$ mass eigenstates, $`\mathrm{\Delta }t`$ = $`t_{CP}`$ $``$ $`t_{\mathrm{tag}}`$, and the $`b`$-flavor charge $`q`$ = +1 ($`1`$) when the tagging $`B`$ meson is a $`B^0`$ ($`\overline{B}^0`$). The $`CP`$-violating parameters $`𝒮`$ and $`𝒜`$ are given by

$$𝒮\frac{2m(\lambda )}{|\lambda |^2+1},𝒜\frac{|\lambda |^21}{|\lambda |^2+1},$$
(2)
where $`\lambda `$ is a complex parameter that depends on both the $`B^0\overline{B}^0`$ mixing and on the amplitudes for $`B^0`$ and $`\overline{B}^0`$ decay to $`f_{CP}`$. To a good approximation in the SM, $`|\lambda |`$ is equal to the absolute value of the ratio of the $`\overline{B}{}_{}{}^{0}f_{CP}`$ to $`B^0f_{CP}`$ decay amplitudes. The SM predicts $`𝒮=\xi _f\mathrm{sin}2\varphi _1`$, where $`\xi _f=+1(1)`$ corresponds to $`CP`$-even (-odd) final states; and $`𝒜=0`$ (or equivalently $`|\lambda |=1`$) for both $`bc\overline{c}s`$ and $`bs\overline{s}s`$ transitions.

In this paper, we report the first measurement of $`CP`$ asymmetries in $`B^0\varphi K_S^0`$ and $`K^+K^{}K_S^0`$ decays, and an improved measurement for $`B^0\eta ^{}K_S^0`$ decay bib:Belle\_etapks  based on a $`78\mathrm{fb}^1`$ data sample, which contains 85 million $`B\overline{B}`$ pairs. Data are collected with the Belle detector at the KEKB energy-asymmetric $`e^+e^{}`$ (3.5 on 8 GeV) collider bib:KEKB  operating at the $`\mathrm{{\rm Y}}(4S)`$ resonance. At KEKB, the $`\mathrm{{\rm Y}}(4S)`$ is produced with a Lorentz boost of $`\beta \gamma =0.425`$ nearly along the electron beamline ($`z`$). Since the $`B^0`$ and $`\overline{B}^0`$ mesons are approximately at rest in the $`\mathrm{{\rm Y}}(4S)`$ center-of-mass system (cms), $`\mathrm{\Delta }t`$ can be determined from the displacement in $`z`$ between the $`f_{CP}`$ and $`f_{\mathrm{tag}}`$ decay vertices: $`\mathrm{\Delta }t(z_{CP}z_{\mathrm{tag}})/\beta \gamma c\mathrm{\Delta }z/\beta \gamma c`$.

The Belle detector bib:Belle  is a large-solid-angle spectrometer that includes a three-layer silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), time-of-flight (TOF) scintillation counters, and an electromagnetic calorimeter comprised of CsI(Tl) crystals (ECL) located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside of the coil is instrumented to detect $`K_L^0`$ mesons and to identify muons (KLM).

We reconstruct $`B^0`$ decays to $`\varphi K_S^0`$ and $`\eta ^{}K_S^0`$ final states for $`\xi _f=1`$, and $`B^0K^+K^{}K_S^0`$ decays that are a mixture of $`\xi _f=+1`$ and $`1`$. $`K^+K^{}`$ pairs that are consistent with $`\varphi K^+K^{}`$ decay are excluded from the $`B^0K^+K^{}K_S^0`$ sample. We find that the $`K^+K^{}K_S^0`$ state is primarily $`\xi _f=+1`$; the $`\xi _f=+1`$ fraction is $`1.04\pm 0.19\text{(stat)}\pm 0.06\text{(syst)}`$ bib:Garmash . In the following determination of $`𝒮`$ and $`𝒜`$, we fix $`\xi _f=+1`$. The intermediate meson states are reconstructed from the following decay chains: $`\eta ^{}\rho ^0(\pi ^+\pi ^{})\gamma `$ or $`\eta ^{}\pi ^+\pi ^{}\eta (\gamma \gamma )`$, $`K_S^0\pi ^+\pi ^{}`$, and $`\varphi K^+K^{}`$.

Candidate $`K_S^0\pi ^+\pi ^{}`$ decays are oppositely charged track pairs that have an invariant mass within 15 (12) MeV/$`c^2`$ of the nominal $`K_S^0`$ mass for the $`B^0\varphi K_S^0`$ ($`B^0K^+K^{}K_S^0`$) mode. The displacement of the $`\pi ^+\pi ^{}`$ vertex from the nominal interaction point (IP) in the plane transverse to the positron beam axis ($`r`$-$`\varphi `$ plane) is required to be greater than 0.1 cm and less than 20 cm. The direction of the combined pion-pair momentum in the $`r`$-$`\varphi `$ plane is required to be within 0.2 radians of the direction defined by the IP and the displaced vertex. We select candidate $`\varphi K^+K^{}`$ decays requiring that the $`K^+K^{}`$ invariant mass is within 10 MeV/$`c^2`$ of the nominal $`\varphi `$ meson mass, the $`\varphi `$ meson momentum in the cms exceeds 2.0 GeV/$`c`$, and the $`K^+K^{}`$ vertex is consistent with the IP. Since the $`\varphi `$ meson selection is effective in reducing background events, we impose only minimal kaon-identification requirements. For selection of non-resonant $`K^+K^{}K_S^0`$ candidates, more stringent kaon-identification requirements, which retain 86% of $`K^\pm `$ at a 7% fake rate for $`\pi ^\pm `$, are used. In this case, charged tracks that are positively identified as electrons or protons are excluded. In addition to the rejection of $`\varphi `$ meson candidates, we reject $`K^+K^{}`$ pairs that are consistent with $`D^0K^+K^{}`$ or $`\chi _{c0}K^+K^{}`$ decay.

For reconstructed $`Bf_{CP}`$ candidates, we identify $`B`$ meson decays using the energy difference $`\mathrm{\Delta }EE_B^{\mathrm{cms}}E_{\mathrm{beam}}^{\mathrm{cms}}`$ and the beam-energy constrained mass $`M_{\mathrm{bc}}\sqrt{(E_{\mathrm{beam}}^{\mathrm{cms}})^2(p_B^{\mathrm{cms}})^2}`$, where $`E_{\mathrm{beam}}^{\mathrm{cms}}`$ is the beam energy in the cms, and $`E_B^{\mathrm{cms}}`$ and $`p_B^{\mathrm{cms}}`$ are the cms energy and momentum of the reconstructed $`B`$ candidate, respectively. The $`B`$ meson signal region is defined as $`5.27\mathrm{GeV}/c^2<M_{\mathrm{bc}}<5.29\mathrm{GeV}/c^2`$ and $`|\mathrm{\Delta }E|<0.051`$ GeV for $`B^0\varphi K_S^0`$ or $`|\mathrm{\Delta }E|<0.040`$ GeV for $`B^0K^+K^{}K_S^0`$. In order to suppress background from the $`e^+e^{}q\overline{q}`$ ($`q=u,d,s`$) and $`c\overline{c}`$ continuum, we form signal and background likelihood functions, $`_\mathrm{S}`$ and $`_{\mathrm{BG}}`$, from a set of variables that characterize the event topology bib:Garmash ; bib:Belle\_etapks . We determine $`_\mathrm{S}`$ from Monte Carlo (MC) and $`_{\mathrm{BG}}`$ from data, and impose mode-dependent thresholds on the likelihood ratio $`_\mathrm{S}/(_\mathrm{S}+_{\mathrm{BG}})`$. The numbers of reconstructed candidates are 59 and 230 for $`B^0\varphi K_S^0`$ and for $`B^0K^+K^{}K_S^0`$, respectively.

For $`B^0\eta ^{}K_S^0`$ decay, we use the same selection criteria as those used in our previously published analysis bib:Belle\_etapks  if both of the charged pions in the $`\eta ^{}\pi ^+\pi ^{}\eta `$ or the $`\eta ^{}\rho ^0\gamma `$ decay have associated SVD hits. We also reconstruct events where only one of the charged pions has associated SVD hits. In this case, the requirement on the impact parameter is relaxed for the track without SVD hits, while a higher threshold is imposed on the likelihood ratio. The number of reconstructed $`B^0\eta ^{}K_S^0`$ candidates is 311.

Charged leptons, kaons, pions, and $`\mathrm{\Lambda }`$ baryons that are not associated with the reconstructed $`f_{CP}`$ decay are used to identify the $`b`$-flavor of the accompanying $`B`$ meson, which decays into $`f_{\mathrm{tag}}`$. Based on the measured properties of these tracks, two parameters, $`q`$ and $`r`$, are assigned to each event. The first, $`q`$, has the discrete value $`+1`$ ($`1`$) when the tag-side $`B`$ meson is more likely to be a $`B^0`$ ($`\overline{B}^0`$). The parameter $`r`$ is an event-by-event MC-determined flavor-tagging dilution factor that ranges from $`r=0`$ for no flavor discrimination to $`r=1`$ for an unambiguous flavor assignment. It is used only to sort data into six intervals of $`r`$, according to the estimated flavor purity. The wrong-tag probabilities for each of these intervals, $`w_l(l=1,6)`$, which are used in the final fit, are determined directly from the data. Samples of $`B^0`$ decays to exclusively reconstructed self-tagging channels are utilized to obtain $`w_l`$ using time-dependent $`B^0`$-$`\overline{B}^0`$ mixing: $`(N_{\mathrm{OF}}N_{\mathrm{SF}})/(N_{\mathrm{OF}}+N_{\mathrm{SF}})=(12w_l)\mathrm{cos}(\mathrm{\Delta }m_d\mathrm{\Delta }t)`$, where $`N_{\mathrm{OF}}`$ and $`N_{\mathrm{SF}}`$ are the numbers of opposite ($`B^0\overline{B}{}_{}{}^{0}B^0\overline{B}^0`$) and same ($`B^0\overline{B}{}_{}{}^{0}B^0B^0,\overline{B}{}_{}{}^{0}\overline{B}^0`$) flavor events. The event fractions and wrong tag fractions for each $`r`$ interval are described elsewhere bib:Belle\_sin2phi1\_78fb-1 .

The decay vertices of $`B^0`$ mesons are reconstructed using tracks that have enough SVD hits: i.e. both $`z`$ and $`r`$-$`\varphi `$ hits in at least one SVD layer and at least one additional layer with a $`z`$ hit, where the $`r`$-$`\varphi `$ plane is perpendicular to the $`z`$ axis. Each vertex position is required to be consistent with the IP profile, which is determined run-by-run and smeared in the $`r`$-$`\varphi `$ plane by 21 $`\mu `$m to account for the $`B`$ meson decay length. With these requirements, we are able to determine a vertex even with a single track. The vertex position for the $`f_{CP}`$ decay is reconstructed using charged kaons for $`B^0\varphi K_S^0`$ and $`B^0K^+K^{}K_S^0`$ decays and using charged pions from $`\rho ^0`$ or $`\eta ^{}`$ decays for $`B^0\eta ^{}K_S^0`$. The algorithm for the $`f_{\mathrm{tag}}`$ vertex reconstruction is chosen to minimize the effect of long-lived particles, secondary vertices from charmed hadrons and a small fraction of poorly reconstructed tracks bib:resol\_nim . From all the charged tracks with associated SVD hits except those used for $`f_{CP}`$, we select tracks with a position error in the $`z`$ direction of less than 500 $`\mu `$m, and with an impact parameter with respect to the $`f_{CP}`$ vertex of less than 500 $`\mu `$m. Track pairs with opposite charges are removed if they form a $`K_S^0`$ candidate with an invariant mass within $`\pm 15`$ MeV/$`c^2`$ of the nominal $`K_S^0`$ mass. If the reduced $`\chi ^2`$ associated with the $`f_{\mathrm{tag}}`$ vertex exceeds 20, the track making the largest $`\chi ^2`$ contribution is removed and the vertex is refitted. This procedure is repeated until an acceptable reduced $`\chi ^2`$ is obtained.

After flavor tagging and vertex reconstruction, we obtain the numbers of $`B^0f_{CP}`$ candidates, $`N_{\mathrm{ev}}`$, listed in Table 1.

Figure 1 shows the $`M_{\mathrm{bc}}`$ distributions for the reconstructed $`B`$ candidates that have $`\mathrm{\Delta }E`$ values within the signal region.

To assign an event-by-event signal probability for use in the maximum-likelihood fit of the $`CP`$-violating parameters, we determine event distribution functions in the $`\mathrm{\Delta }E`$-$`M_{\mathrm{bc}}`$ plane for both signal and background. The signal distribution is modeled with a single two-dimensional Gaussian, where the widths are allowed to float in the fit to data. For the continuum background, we use a linear function for $`\mathrm{\Delta }E`$ and the ARGUS parameterization bib:ARGUS  for $`M_{\mathrm{bc}}`$. We use events outside the signal region as well as a large MC sample to study the background components. The dominant background comes from continuum events. In addition, according to MC simulation, there is a non-negligible ($`8\%`$) contamination from $`B\overline{B}`$ background events in $`B^0\eta ^{}K_S^0(\eta ^{}\rho ^0\gamma )`$. The contributions from $`B\overline{B}`$ events are smaller for other decay modes. The contamination of $`K^+K^{}K_S^0`$ events in the $`\varphi K_S^0`$ sample (and vice versa) is also small and is treated as a source of systematic uncertainty. Finally, backgrounds from $`B^0f_0(980)K_S^0`$ decay, which has the opposite $`CP`$ eigenvalue to $`\varphi K_S^0`$, are found to be negligible.

We determine $`𝒮`$ and $`𝒜`$ for each mode by performing an unbinned maximum-likelihood fit to the observed $`\mathrm{\Delta }t`$ distribution. The probability density function (PDF) expected for the signal distribution is given by Eq. (1) with $`q`$ replaced by $`q(12w_l)`$ to account for the effect of incorrect flavor assignment. The distribution is convolved with the proper-time interval resolution function $`R_{\mathrm{sig}}(\mathrm{\Delta }t)`$, which takes into account the finite vertex resolution. It is formed by convolving four components: the detector resolutions for $`z_{CP}`$ and $`z_{\mathrm{tag}}`$, the shift in the $`z_{\mathrm{tag}}`$ vertex position due to secondary tracks originating from charmed particle decays, and the kinematic approximation that the $`B`$ mesons are at rest in the cms bib:resol\_nim . A small component of broad outliers in the $`\mathrm{\Delta }z`$ distribution, caused by mis-reconstruction, is represented by a Gaussian function $`P_{\mathrm{ol}}(\mathrm{\Delta }t)`$. We determine twelve resolution parameters and the neutral- and charged-$`B`$ lifetimes simultaneously from a fit to the $`\mathrm{\Delta }t`$ distributions of hadronic $`B`$ decays and obtain an average $`\mathrm{\Delta }t`$ resolution of $`1.43`$ps (rms). We determine the following likelihood value for each event:

$`P_i(\mathrm{\Delta }t_i;𝒮,𝒜)=(1f_{\mathrm{ol}}){\displaystyle _{\mathrm{}}^{\mathrm{}}}[f_{\mathrm{sig}}𝒫_{\mathrm{sig}}(\mathrm{\Delta }t^{},q,w_l)R_{\mathrm{sig}}(\mathrm{\Delta }t_i\mathrm{\Delta }t^{})`$
$`+(1f_{\mathrm{sig}})𝒫_{\mathrm{bkg}}(\mathrm{\Delta }t^{})R_{\mathrm{bkg}}(\mathrm{\Delta }t_i\mathrm{\Delta }t^{})]d(\mathrm{\Delta }t^{})+f_{\mathrm{ol}}P_{\mathrm{ol}}(\mathrm{\Delta }t_i)`$ (3)
where $`f_{\mathrm{ol}}`$ is the outlier fraction and $`f_{\mathrm{sig}}`$ is the signal probability calculated as a function of $`\mathrm{\Delta }E`$ and $`M_{\mathrm{bc}}`$. $`𝒫_{\mathrm{bkg}}(\mathrm{\Delta }t)`$ is a PDF for background events, which dilutes the significance of $`CP`$ violation in Eq. (1). It is modeled as a sum of exponential and prompt components, and is convolved with a sum of two Gaussians, $`R_{\mathrm{bkg}}`$, which represents a resolution function for the background. All parameters in $`𝒫_{\mathrm{bkg}}(\mathrm{\Delta }t)`$ and $`R_{\mathrm{bkg}}`$ are determined by the fit to the $`\mathrm{\Delta }t`$ distribution of a background-enhanced control sample bib:BBbg ; i.e. events away from the $`\mathrm{\Delta }E`$-$`M_{\mathrm{bc}}`$ signal region. We fix the $`\tau _{B^0}`$ and $`\mathrm{\Delta }m_d`$ at their world-average values bib:PDG . The only free parameters in the final fit are $`𝒮`$ and $`𝒜`$, which are determined by maximizing the likelihood function

$$L=\underset{i}{}P_i(\mathrm{\Delta }t_i;𝒮,𝒜)$$
(4)
where the product is over all events. Table 2 summarizes the results of the fit. The table shows the values of $`𝒜`$ and $`\xi _f𝒮`$, which, in the SM, is equal to $`\mathrm{sin}2\varphi _1`$.

The first errors are statistical and the second errors are systematic. The third error for the $`K^+K^{}K_S^0`$ mode arises from the uncertainty in the fraction of the $`CP`$-odd component bib:Garmash . Figure 2 shows the observed $`\mathrm{\Delta }t`$ distribution for $`q\xi _f=1`$ (upper figure) and $`q\xi _f=+1`$ (lower figure) event samples for each $`f_{CP}`$ mode.

Figure 3 shows the raw asymmetry in each $`\mathrm{\Delta }t`$ bin without background subtraction, which is defined by

$$A\frac{N_{q\xi _f=1}N_{q\xi _f=+1}}{N_{q\xi _f=1}+N_{q\xi _f=+1}},$$
(5)
where $`N_{q\xi _f=+1(1)}`$ is the number of observed candidates with $`q\xi _f=+1(1)`$. The curves show the results of the unbinned-maximum likelihood fit to the asymmetry distribution, $`\xi _f𝒮\mathrm{sin}(\mathrm{\Delta }m_d\mathrm{\Delta }t)\xi _f𝒜\mathrm{cos}(\mathrm{\Delta }m_d\mathrm{\Delta }t)`$. These are the first measurements of the $`CP`$ violation parameters for $`B^0\varphi K_S^0`$ and $`B^0K^+K^{}K_S^0`$ decays. The result for $`\eta ^{}K_S^0`$ supersedes the previous result bib:Belle\_etapks . We obtain values consistent with the present world average of $`\mathrm{sin}2\varphi _1=+0.734\pm 0.054`$ bib:Nir02  in $`B^0K^+K^{}K_S^0`$ and $`\eta ^{}K_S^0`$ decays, while a 2.1$`\sigma `$ deviation is observed in $`B^0\varphi K_S^0`$ decay.

Fits to the same samples with the direct $`CP`$ violation parameter $`𝒜`$ fixed at zero yield $`\xi _f𝒮=0.83\pm 0.72`$(stat) for $`B^0\varphi K_S^0`$, $`\xi _f𝒮=+0.59\pm 0.47`$(stat) for $`B^0K^+K^{}K_S^0`$, and $`\xi _f𝒮=+0.77\pm 0.38`$(stat) for $`B^0\eta ^{}K_S^0`$. As a consistency check for the $`𝒮`$ term, we select the charged $`B`$ meson decays $`B^+\varphi K^+`$ and $`B^+\eta ^{}K^+`$ and apply the same fit procedure. We obtain $`𝒮=0.05\pm 0.32`$(stat), $`𝒜=0.29\pm 0.24`$(stat) for $`B^+\varphi K^+`$ decay and $`𝒮=0.03\pm 0.20`$(stat), $`𝒜=0.05\pm 0.13`$(stat) for $`B^+\eta ^{}K^+`$ decay. Both results on the $`𝒮`$ term are consistent with no $`CP`$ asymmetry, as expected.

The largest source of systematic error for the $`B^0\varphi K_S^0`$ mode is the uncertainty in the signal fraction and the background $`\mathrm{\Delta }t`$ shape ($`\pm 0.17`$ for $`𝒮`$ and $`\pm 0.14`$ for $`𝒜`$ in total) determined from the events in the sideband regions in the $`\mathrm{\Delta }E`$-$`M_{\mathrm{bc}}`$ plane. Other significant contributions come from uncertainties in the vertex reconstruction, the resolution function parameters, wrong tag fractions, $`\tau _{B^0}`$, and $`\mathrm{\Delta }m_d`$. We add each contribution in quadrature to obtain the total systematic uncertainty. Systematic uncertainties from these sources are also examined for the other modes. We find that the largest uncertainties arise from the vertex reconstruction ($`\pm 0.09`$ for $`𝒮`$ and $`\pm 0.08`$ for $`𝒜`$) for the $`B^0K^+K^{}K_S^0`$ mode, and from the resolution function parameters ($`{}_{0.04}{}^{}{}_{}{}^{+0.03}`$ for $`𝒮`$) and the signal fraction ($`\pm 0.02`$ for $`𝒜`$) for the $`B^0\eta ^{}K_S^0`$ mode.

In summary, we have performed the first measurement of $`CP`$ violation parameters in the $`B^0\varphi K_S^0`$ and $`K^+K^{}K_S^0`$ decays. We also provide an improved measurement for $`\eta ^{}K_S^0`$ decay. These modes are dominated by the $`bs\overline{s}s`$ transition and are sensitive to possible new $`CP`$-violating phases. Our results for $`B^0\eta ^{}K_S^0`$ and $`K^+K^{}K_S^0`$ are consistent with those obtained for $`B^0J/\psi K_S^0`$ and other decays governed by the $`bc\overline{c}s`$ transition. A $`2.1\sigma `$ deviation is observed for $`B^0\varphi K_S^0`$.

###### Acknowledgements.

We wish to thank the KEKB accelerator group for the excellent operation of the KEKB accelerator. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology of Japan and the Japan Society for the Promotion of Science; the Australian Research Council and the Australian Department of Industry, Science and Resources; the National Science Foundation of China under contract No. 10175071; the Department of Science and Technology of India; the BK21 program of the Ministry of Education of Korea and the CHEP SRC program of the Korea Science and Engineering Foundation; the Polish State Committee for Scientific Research under contract No. 2P03B 17017; the Ministry of Science and Technology of the Russian Federation; the Ministry of Education, Science and Sport of the Republic of Slovenia; the National Science Council and the Ministry of Education of Taiwan; and the U.S. Department of Energy.


format: text

text:
# Spectral pollution and second order relative spectra for self-adjoint operators

## 1 Introduction

Let $`A`$ be a self-adjoint operator on a Hilbert space $``$. In order to determine the spectrum $`\mathrm{Spec}(A)`$ of $`A`$, one can use the following projection method. Let $`(_k)_k`$ be a sequence of closed linear subspaces of the domain $`\mathrm{Dom}(A)`$ of $`A`$ and suppose that the corresponding orthogonal projections $`P_k:_k`$ converge strongly to the identity operator $`I`$. Let $`\mathrm{\Lambda }(A)`$ be the set of all such sequences of subspaces. One may ask whether or not

$$\underset{k\mathrm{}}{lim}\mathrm{Spec}(A,_k)=\mathrm{Spec}(A),$$
(1)
where $`\mathrm{Spec}(A,_k)`$ is the spectrum of $`P_kA:_k_k`$ and “$`lim`$” is defined in an appropriate sense. Unfortunately the answer is, in general, negative and $`\underset{k\mathrm{}}{lim}\mathrm{Spec}(A,_k)`$ may contain points which do not belong to $`\mathrm{Spec}(A)`$. This phenomenon is well known in numerical analysis and is called spectral pollution.

Since we allow $`A`$ to be unbounded, it is convenient for us to extend $`\mathrm{Spec}(A)`$ by adding to it $`\mathrm{}`$, $`+\mathrm{}`$, or both, if $`A`$ is unbounded below, above or from both sides, respectively. We denote the resulting set by $`\widehat{\mathrm{S}}\mathrm{pec}(A)`$. In other words, $`\widehat{\mathrm{S}}\mathrm{pec}(A)`$ is the closure of $`\mathrm{Spec}(A)`$ in the two point compactification $`\overline{}:=\{\pm \mathrm{}\}`$ of $``$. Let $`\mathrm{Spec}_{\text{ess}}(A)`$ be the essential spectrum of $`A`$, i.e.

$$\mathrm{Spec}_{\text{ess}}(A):=\mathrm{Spec}(A)\{\text{isolated eigenvalues of finite multiplicity}\},$$

and let $`\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A)`$ be the extended essential spectrum of $`A`$, i.e.

$$\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A):=\widehat{\mathrm{S}}\mathrm{pec}(A)\{\text{isolated eigenvalues of finite multiplicity}\}.$$

It turns out that spectral pollution may occur at any point lying in a gap of the extended essential spectrum, i.e. at any point of the set $`\text{conv}\left(\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A)\right)\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A)`$, where “conv” means the convex hull (see the proof of Theorem 2.1 below, where this fact is demonstrated by choosing an appropriate sequence of finite-dimensional subspaces).

One can argue that the sequence $`(_k)_k`$ in the proof of Theorem 2.1 had been deliberately chosen in a particularly nasty way and that things like this do not usually happen in “real life”. This argument is only partially true at best, and examples in Sections 3 and 4 show how the standard approach leads to spectral pollution in some very straightforward problems.

One can try to deal with spectral pollution by choosing $`_k`$ in a “reasonable way”. However the known recipes (see, e.g., \[BBG\], \[BDG\], \[RSHSPV\], and references therein) are by no means universal or rigorously justified for a wide range of problems. An absolutely “safe” choice of $`(_k)_k`$ would be the one based on the spectral decomposition of $`A`$. Unfortunately knowing the spectral decomposition is much more than knowing just the spectrum which is to be found.

We propose a complementary (not an alternative!) way of dealing with spectral pollution, which is based on the notion of the second order spectrum $`\mathrm{Spec}_2(A,)`$ of $`A`$ relative to a subspace $`\mathrm{Dom}(A)`$ (see Definition 2.3) introduced in \[Dav\]. Second order relative spectra usually contain complex numbers, but it turns out (see Theorem 2.5) that if $`z\mathrm{Spec}_2(A,)`$ then

$$\mathrm{Spec}(A)[\mathrm{Re}z|\mathrm{Im}z|,\mathrm{Re}z+|\mathrm{Im}z|]\mathrm{}.$$
(2)
In particular, this means that if $`z`$ belongs to a second order relative spectrum of $`A`$, and its imaginary part is small, than there is a point of the spectrum of $`A`$ close to $`\mathrm{Re}z`$.

This fact suggests the following strategy, which uses (2) as a posteriori estimate. One should choose a reasonable sequence $`(_k)_k`$ and apply the standard projection method, i.e. find $`\mathrm{Spec}(A,_k)`$. Additionally, one should find $`\mathrm{Spec}_2(A,\stackrel{~}{}_k)`$ (for possibly different $`\stackrel{~}{}_k\mathrm{Dom}(A)`$). Suppose $`\mathrm{Spec}(A,_k)`$ contain points which appear to converge to some $`\lambda `$ as $`k\mathrm{}`$. If $`\mathrm{Spec}_2(A,\stackrel{~}{}_k)`$ also contain points lying close to $`\lambda `$, then $`\lambda `$ is indeed close to the spectrum of $`A`$. If however $`\mathrm{Spec}_2(A,\stackrel{~}{}_k)`$ do not contain points approximating $`\lambda `$, one cannot conclude that $`\lambda `$ is far from the spectrum. A definite claim which can be made in this case is that $`\lambda `$ requires additional attention111One can, for example, use the method suggested in the companion paper by E.B. Davies and M. Plum \[DP\]., as spectral pollution may be around the corner. It has to be said though, that our (very limited!) experience seems to suggest that $`\mathrm{Spec}_2(A,\stackrel{~}{}_k)`$ do capture the spectrum quite well.

We do not claim that the above method is going to be universally successful. Our main ambition is merely to bring it to the attention of those concerned with finding spectra numerically. We have applied it to two problems which were chosen mainly for their simplicity. The results (which turned out to be better than we expected) are presented in Sections 3, 4, 6, and 7.

Section 2 contains the main theoretical results on standard (first order) and second order relative spectra. Similar results on second (and higher) order relative spectra have been obtained in \[Sha\] under slightly stronger restrictions. We start with Theorem 2.1, which has been mentioned above and which says that spectral pollution may occur at any point lying in a gap of the extended essential spectrum. It implies that one should be particularly careful when applying a projection method to find the spectrum of an operator $`A`$ in the following cases:
a) $`A`$ is not semi-bounded and $`\mathrm{Spec}(A)`$;
b) $`A`$ is semi-bounded but unbounded, and has a nonempty essential spectrum different from a half-line;
c) $`A`$ is bounded and has at least one gap in the essential spectrum.
Examples of types c) and b) are treated in Sections 3 and 4 respectively. We return to these examples in Sections 6 and 7 after giving the proofs of the main theoretical results in Section 5.

Section 2 also discusses an additional difficulty which arises for unbounded operators. It is well known that if $`A`$ is bounded, then the left-hand side of (1) contains the right-hand side. If, however, $`A`$ is unbounded above, say, then the left-hand side of (1) may well be equal to $`\{+\mathrm{}\}`$ (see Theorem 2.2).

In order to make the paper more readable, we have divided the material into two parts. The first one consists of Sections 2–4 and covers the main theoretical and numerical results. The same material is treated at a greater depth in the second part of the paper, which consists of Sections 5–7. The summary is given in Section 8.

Finally, the accompanying webpage

www.ma.hw.ac.uk/$``$levitin/SpectralPollution

contains the listings of the Matlab programmes written by the authors and used for numerical calculations in Examples I and II.

## 2 Main results

Let $`A`$ be a self-adjoint operator on a Hilbert space $``$. Throughout the paper $`E()`$ denote the spectral projections associated with $`A`$. The proofs of all results of this Section are given in Section 5.

The following result states that spectral pollution may occur at any point in a gap of the extended essential spectrum of $`A`$.

###### Theorem 2.1.

For any $`\lambda \text{conv}\left(\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A)\right)\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A)`$ there exists an increasing sequence $`(_k)_k\mathrm{\Lambda }(A)`$ such that

$$\lambda \mathrm{Spec}(A,_k),k.$$

It is well known that for bounded self-adjoint operators the left-hand side of (1) contains the right-hand side:

$$\underset{k\mathrm{}}{lim}\mathrm{Spec}(A,_k)\mathrm{Spec}(A)$$
(3)
(see, e.g., \[Arv\] or \[Sha\]). So, spectral pollution means that the left-hand side of (1) is strictly larger than the right-hand side. For unbounded self-adjoint operators the situation is more complicated: the left-hand side of (1) does not necessarily contain the right-hand side. The reason for this is intuitively very clear. Let $`_k`$ be finite dimensional and suppose for definiteness that $`A`$ is not bounded above. Then perturbing a basis of $`_k`$ by arbitrarily small vectors lying in $`\mathrm{Dom}(A)E([N,+\mathrm{}))`$ with sufficiently large $`N`$, one can push the spectrum of the corresponding finite dimensional operator arbitrarily far to the right.

###### Theorem 2.2.

Let $`A`$ be an unbounded above self-adjoint operator on a separable Hilbert space $``$. Then for an arbitrary sequence $`(_k)_k\mathrm{\Lambda }(A)`$ of finite dimensional subspaces and arbitrary sequences $`(\epsilon _k)_k0`$ and $`(R_k)_k+\mathrm{}`$, there exists a sequence $`(_k^{})_k\mathrm{\Lambda }(A)`$ such that

$$P_kP_k^{}<\epsilon _k\text{and}\mathrm{Spec}(A,_k^{})(R_k,+\mathrm{}),k,$$

where $`P_k^{}:_k^{}`$ are the corresponding orthogonal projections. A similar statement holds for operators unbounded below.

The requirement $`_k\mathrm{Dom}(A)`$ is too restrictive for many applications, in particular for approximating spectra of differential operators by finite element methods. When using finite element methods, one usually states the problem in a weak (variational) form and requires that $`_k`$ belong to the domain of the corresponding sesquilinear form. The restriction $`_k\mathrm{Dom}(A)`$ would mean extra smoothness of the finite elements.

Suppose $`A`$ is a bounded below self-adjoint operator on a separable Hilbert space $``$ and let $`Q_A:\times `$ be the corresponding closed sesquilinear form (see, e.g., \[Kat, §1, Ch. VI\]). Let $`(_k)_k`$ be a sequence of finite dimensional linear subspaces of the domain $`\mathrm{Dom}(Q_A)`$ of $`Q_A`$ such that the corresponding orthogonal projections $`P_k:_k`$ converge strongly to $`I`$. Let $`\mathrm{\Lambda }(Q_A)`$ be the set of all such sequences. A number $`\lambda `$ is said to belong to $`\mathrm{Spec}(A,_k)`$ if there exists $`u_k\{0\}`$ such that

$$Q_A(u,v)\lambda (u,v)=0,v_k.$$

If $`_k\mathrm{Dom}(A)`$, then $`Q_A(u,v)=(Au,v)`$ and the above definition means that
$``$ there exists $`u_k\{0\}`$ such that $`(Au,v)\lambda (u,v)=0`$, $`v_k`$,
$``$ i.e. there exists $`u_k\{0\}`$ such that $`P_kAu=\lambda u`$
$``$ i.e. $`\lambda `$ belongs to the spectrum of $`P_kA:_k_k`$.
Hence one gets the old definition of $`\mathrm{Spec}(A,_k)`$ for finite dimensional $`_k\mathrm{Dom}(A)`$ (cf. Remarks 2.4 and 2.6(b) below).

It is not difficult to prove that Theorem 2.2 remains true if $`\mathrm{\Lambda }(A)`$ is replaced by $`\mathrm{\Lambda }(Q_A)`$ (see Section 5).

The above results show the difficulties one may face when trying to approximate the spectrum of a self-adjoint operator with the help of a projection method. A possible way of addressing spectral pollution is based on the notion of the second order relative spectrum.

###### Definition 2.3.

Let $`A`$ be a self-adjoint operator on a Hilbert space $``$ and let $``$ be a finite dimensional subspace of $`\mathrm{Dom}(A)`$. A number $`z`$ is said to belong to the second order spectrum $`\mathrm{Spec}_2(A,)`$ of $`A`$ relative to $``$ if there exists $`u\{0\}`$ such that

$$((AzI)u,(A\overline{z}I)v)=0,v,$$
(4)
where $`(,)`$ is the scalar product in $``$.

###### Remark 2.4.

Suppose $`\mathrm{Dom}(A^2)`$ and let $`P`$ be the orthogonal projection onto $``$. Then $`z\text{Spec}_2(A,)`$
$``$ iff there exists $`u\{0\}`$ such that $`((AzI)^2u,v)=0`$, $`v`$,
$``$ i.e. iff there exists $`u\{0\}`$ such that $`u\text{Ker}(P(AzI)^2)`$,
$``$ i.e. iff the operator $`P(AzI)^2:`$ is not invertible.
The last is the definition of the second order relative spectrum introduced in \[Dav\] and used in \[Sha\].

###### Theorem 2.5.

If $`z\mathrm{Spec}_2(A,)`$ then

$$\mathrm{Spec}(A)[\mathrm{Re}z|\mathrm{Im}z|,\mathrm{Re}z+|\mathrm{Im}z|]\mathrm{}.$$

###### Remark 2.6.

(a) Theorem 2.5 has been proved in \[Sha\] in the case $`\mathrm{Dom}(A^2)`$.
(b) Everything said above about second order relative spectra remains unchanged if $``$ is infinite dimensional, with only one exception: the implication

$$P(AzI)^2:\text{ is not invertible }\text{Ker}(P(AzI)^2)\{0\}$$

(see Remark 2.4) is no longer true.

## 3 Example I: operator of multiplication by a discontinuous function

Consider the operator of multiplication by a function $`a`$:

$$A=aI:L_2([\pi ,\pi ])L_2([\pi ,\pi ]),$$

where

$$a(x)=\{\begin{array}{cc}\frac{3}{2}+\frac{1}{2}\mathrm{cos}\sqrt{5}x,\hfill & \text{for}\pi x<0,\hfill \\ 2+\mathrm{cos}\sqrt{2}x,\hfill & \text{for}0x<\pi .\hfill \end{array}$$

It is clear that

$$\mathrm{Spec}(A)=\mathrm{Spec}_{\text{ess}}(A)=a\left([\pi ,\pi ]\right)=[2,1][1,3].$$

The function $`a`$ and, in particular, the irrational factors $`\sqrt{5}`$ and $`\sqrt{2}`$ have been chosen to get a “generic” situation and avoid any possible hidden symmetries.

Using the standard orthonormal basis $`e_k(x):=(2\pi )^{1/2}\mathrm{exp}(ikx)`$, $`k`$, one can represent $`A`$ as a (doubly) infinite Toeplitz matrix and take the orthogonal projections onto $`_N:=\text{span}(\{e_k\}_{k=n}^n)`$, where $`N=\text{dim}_N=2n+1`$. Figure 2 shows the approximation of the spectrum obtained by the standard projection method for a typical value of $`N`$. One can see that the essential spectrum $`[2,1][1,3]`$ is approximated very well, but that, additionally, there are some computed eigenvalues lying in the gap of the essential spectrum.

The spurious eigenvalues vary very slowly with $`N`$, see Table 1, and their number $`\mathrm{\#}(\mathrm{Spec}(A,_N)(1,1))`$ also increases slowly as $`N`$ increases. In fact, they seem to converge monotonically to some points in $`(1,1)`$, and judging by the numerical evidence alone, one might have thought that the limits belonged to the spectrum, if the example was less trivial and the spectrum was not known.

The presence of spurious points is hardly surprising. Indeed, if instead of the operator of multiplication by $`a`$, one considers the corresponding Toeplitz operator $`T_a`$ acting on the Hardy space $`𝙷^2=\text{span}(\{e_k\}_{k0})`$, then the projection method using $`_N:=\text{span}(\{e_k\}_{k=0}^{N1})`$ leads to the *same* finite-dimensional Toeplitz matrices. In other words, the standard projection method does not know which spectrum is being calculated: that of the operator of multiplication or of the Toeplitz operator. The latter equals $`[2,3]`$ (see \[HW\] or \[BS, 2.36\]). Therefore not only will we have spectral pollution, the spurious eigenvalues will eventually fill the entire gap $`(1,1)`$ (cf. (3) or \[Böt\]).

###### Remark 3.1.

It is interesting that spurious eigenvalues appear so slowly. Three spurious eigenvalues are present already for rather small values of $`N`$. The fourth one emerges from $`1`$ for $`N161`$, and the fifth one emerges from $`1`$ for $`N631`$. Such a slow rate seems to have the following explanation. The part $`(1,1)`$ of the spectrum of the Toeplitz operator $`T_a`$ is there because of the discontinuities of $`a(x)`$ at $`x=0`$ and $`x=\pi \pi `$. Functions responsible for this part of the spectrum appear to be very localized at $`0`$ and $`\pi \pi `$ and hence have a very slow decay of the Fourier coefficients. Figure 3 shows the graph of the trigonometric polynomial which is an eigenfunction of $`P_NA:_N_N`$, $`N=1001`$, corresponding to $`\lambda 0.4717`$.

We now show how second order relative spectra locate the spectrum of $`A`$ and, in particular, detect spectral pollution in the standard projection method.

Figure 4 shows the second order relative spectra $`\mathrm{Spec}_2(A,_N)`$ computed for $`N=101,201,501`$.

For each $`N`$, the second order relative spectrum consists of a few points lying approximately on a circle centred on the real axis and two contours “surrounding” two intervals of the (essential) spectrum. The number of points of $`\mathrm{Spec}_2(A,_N)`$ lying on the circle does not seem to increase with $`N`$.

Figure 5 gives a detailed view of Figure 4 near the real axis (note the different scales along the coordinate axes). The contours made of the points of the second order relative spectra, which “surround” two intervals of the (essential) spectrum, get closer to the real line as $`N`$ increases. Most importantly, the absence of points of the second order relative spectra close to the eigenvalues found by the standard projection method inside the interval $`(1,1)`$ suggests that these eigenvalues are spurious.

In order to illustrate a situation when an operator has both non-empty essential and discrete spectra, let us consider a two-dimensional perturbation of $`A=aI`$,

$$\begin{array}{cc}\hfill Bu(x):=a(x)u(x)& \frac{3\mathrm{exp}(ix)}{2\pi }_\pi ^\pi u(t)\mathrm{exp}(it)𝑑t\hfill \\ & +\frac{\mathrm{exp}(2ix)}{2\pi }_\pi ^\pi u(t)\mathrm{exp}(2it)𝑑t,x[\pi ,\pi ].\hfill \end{array}$$
(5)
This time, the essential spectrum is again very well approximated by the usual projection method, and the perturbation gives rise to several eigenvalues, see Figure 6, listed here for $`N=401`$ : one ($`\lambda _13.6057`$) to the left of the essential spectrum, six ($`\lambda _20.9842`$, $`\lambda _30.8947`$, $`\lambda _40.8203`$, $`\lambda _50.3496`$, $`\lambda _60.4603`$, and $`\lambda _70.6082`$) in the gap $`(1,1)`$, and one $`\lambda _83.0432`$) to the right of the essential spectrum.

To establish which of the eight eigenvalues are spurious, we compute the second order relative spectra of $`B`$. Their structure is similar to that of $`\mathrm{Spec}_2(A,_N)`$ apart from several additional points located very close to the real axis, see Figures 7 and 8. In particular, for $`N=401`$ these points are approximately $`3.6056\pm 0.0505i`$, $`0.8929\pm 0.0266i`$, $`0.6085\pm 0.0442i`$, and $`3.0413\pm 0.0178i`$. Thus, the eigenvalues $`\lambda _1`$, $`\lambda _3`$, $`\lambda _7`$ and $`\lambda _8`$ are correct, whereas the other eigenvalues $`\lambda _2`$, $`\lambda _4`$, $`\lambda _5`$, and $`\lambda _6`$ are almost certainly spurious.

## 4 Example II: Stokes type system of ODEs

Our next example is an operator $`H`$ in $`L_2([\pi ,\pi ])`$ generated by a simple $`2\times 2`$ Stokes type system of ODEs subject to the periodic boundary conditions,

$$H=\left(\begin{array}{cc}\frac{d^2}{dx^2}& \frac{d}{dx}\\ \frac{d}{dx}& 2\end{array}\right).$$

Finding the spectrum of $`H`$ explicitly is easy. The essential spectrum is found by using a general result on Agmon-Douglis-Nirenberg (that is, mixed order) elliptic systems \[GG\], which gives

$$\mathrm{Spec}_{\text{ess}}(H)=\left\{\lambda \right|\xi \{0\}:det\left(\begin{array}{cc}\xi ^2& i\xi \\ i\xi & 2\lambda \end{array}\right)=0\}=\{1\}.$$

To find the discrete spectrum, consider the system of equations

$$\{\begin{array}{cccc}\hfill u^{\prime \prime }& v^{}\hfill & & =\lambda u,\hfill \\ \hfill u^{}& +2v\hfill & & =\lambda v\hfill \end{array}$$
(6)
subject to the periodic boundary conditions. Expressing $`v`$ from the second equation and substituting the result into the first one gives

$$u^{\prime \prime }=\frac{\lambda (\lambda 2)}{\lambda 1}u.$$
(7)
Comparing this equation with the standard periodic spectral problem for $`d^2/dx^2`$ we conclude that

$$\lambda \mathrm{Spec}_{\text{disc}}(H)\text{iff}\frac{\lambda (\lambda 2)}{\lambda 1}\mathrm{Spec}\left(\frac{d^2}{dx^2}\right)=\{k^2:k\}.$$
(8)
This implies that $`H`$ has two series of eigenvalues $`\lambda _k^\pm `$ given by

$$\lambda _k^\pm =\frac{k^2+2\pm \sqrt{(k^2+2)^24k^2}}{2},k=0,1,2,\mathrm{}$$

Each eigenvalue $`\lambda _k^{}`$ or $`\lambda _k^+`$ with $`k>0`$ is of multiplicity 2. The eigenvalues of the “minus” series are located in the interval $`[0,1)`$ and converge to $`1`$ as $`k\mathrm{}`$. The eigenvalues of the “plus” series are located in the interval $`[2,+\mathrm{})`$ and tend to $`+\mathrm{}`$ as $`k\mathrm{}`$. Note that the interval $`(1,2)`$ is spectrum-free. The eigenvectors $`(u,v)`$ of $`H`$ corresponding to the eigenvalue $`\lambda _k^\pm `$, $`k>0`$ are linear combinations of

$$\left(\begin{array}{c}1\\ ik/(\lambda _k^\pm 2)\end{array}\right)e^{ikx}\text{and}\left(\begin{array}{c}1\\ ik/(\lambda _k^\pm 2)\end{array}\right)e^{ikx}.$$
(9)
The eigenvectors corresponding to $`\lambda _0^{}=0`$ and $`\lambda _0^+=2`$ are proportional to $`\left(\begin{array}{c}1\\ 0\end{array}\right)`$ and $`\left(\begin{array}{c}0\\ 1\end{array}\right)`$ respectively.

The projection method we use in the numerical analysis of the spectrum of $`H`$ is the standard finite element method (FEM). It turns out (see Section 7 below) that the results of calculations depend dramatically upon the choice of the corresponding spaces of finite elements.

We start with the simplest possible choice — the space $`𝐋_N^p`$ of one dimensional periodic Lagrange elements of order $`p`$, that is, of continuous functions which are piecewise polynomials of degree $`p`$. Here $`N`$ is the number of mesh node points on the interval $`[\pi ,\pi ]`$, so the total dimension of $`𝐋_N^p`$ is $`Np`$.

###### Remark 4.1.

In most of the numerical experiments, we use just the linear Lagrange elements (i.e., of order $`1`$). We shall therefore write $`𝐋_N:=𝐋_N^1`$. Additional tests have shown that varying the order affects neither the qualitative picture, nor the rate of convergence.

It is important to note that the spaces $`𝐋_N^p`$ are embedded into the Sobolev space $`W_2^1([\pi ,\pi ])`$ but not into the Sobolev spaces $`W_2^2([\pi ,\pi ])`$. Thus, the elements of $`𝐋_N\times 𝐋_N`$ do not necessarily belong to $`\mathrm{Dom}(H)`$. However, they belong to the domain of the corresponding quadratic form, which is sufficient for the FEM.

The Lagrange elements are provided by all standard FEM packages, including, e.g., FEMLAB (see \[FEM\]).

The calculation of $`\mathrm{Spec}(H,𝐋_N\times 𝐋_N)`$ produces a very striking result: it gives very accurately several initial eigenvalues from the $`\lambda ^\pm `$ series, filling at the same time the entire interval $`(1,2)`$ with spurious eigenvalues (see Figure 9, which, for brevity, shows only the eigenvalues lying in the interval $`[0,6]`$). In the particular case $`N=100`$, there are 15 eigenvalues in $`[0,1]`$, 99 eigenvalues in $`[2,\mathrm{})`$, and the remaining 86 eigenvalues are in $`(1,2)`$.

The subspace $`=𝐋_N\times 𝐋_N`$ used above is unsuitable for the calculation of the second order relative spectra $`\mathrm{Spec}_2(H,)`$ since for elements of this space the left-hand side of (4) does not make sense. Therefore, as usual when dealing with fourth order ODES, let us introduce the standard space of periodic cubic Hermite elements $`𝐌_N`$, whose functions are continuously differentiable and are piecewise polynomials of degree 3. They belong to the Sobolev space $`W_2^2([\pi ,\pi ])`$ so they are smoother than elements of $`𝐋_N`$. Here again $`N`$ is the number of mesh node points, so that $`dim𝐌_N=2N`$.

Since the operator $`H`$ is of mixed order, it is sufficient to choose the space $`𝐌_N\times 𝐋_N`$ for the projection space in (4). The results of calculation of $`\mathrm{Spec}_2(H,𝐌_N\times 𝐋_N)`$ are shown in Figures 10 and 11. The first figure (with different scales along the real and imaginary axes) demonstrates once more the advantage of this method in eliminating the spectral pollution. Note also the series of points of $`\mathrm{Spec}_2(H,𝐌_N\times 𝐋_N)`$ along the line $`\mathrm{Re}\lambda =1`$ — this indicates the presence of the essential spectrum of $`H`$ at $`\lambda =1`$. The second figure shows all computed points of the second order relative spectra for $`N=25,50`$, and $`100`$.

The points of the second order relative spectra give good approximations for the low eigenvalues of $`H`$. We postpone the discussion of the rates of convergence in various cases until Section 7. We note only that the number of well approximated eigenvalues increases with $`N`$, see Figure 11.

Examples I and II are discussed further in Sections 6 and 7; we proceed now to the proofs of the main theoretical results stated in Section 2.

## 5 Proofs of the main results

###### Proof of Theorem 2.1.

It is clear that there exists $`\lambda _{}\widehat{\mathrm{S}}\mathrm{pec}_{\text{ess}}(A)`$ which lies strictly to the left of $`\lambda `$. Consider first the case $`\lambda _{}\mathrm{Spec}_{\text{ess}}(A)`$. If $`\lambda _{}`$ is an isolated point of $`\mathrm{Spec}(A)`$, then it is an eigenvalue of infinite multiplicity and we define $`_j^{}:=\text{span}\left(\{e_j\}\right)`$, where $`e_j`$, $`j`$ are mutually orthogonal eigenvectors of $`A`$ corresponding to $`\lambda _{}`$. If $`\lambda _{}`$ is not an isolated point, there exist pair-wise disjoint compact sets $`W_j^{}`$ lying strictly to the left of $`\lambda `$ and such that $`W_j^{}\mathrm{Spec}(A)\mathrm{}`$, $`j`$. In this case we choose $`_j^{}=E(W_j^{})`$, where $`E()`$ are the spectral projections corresponding to $`A`$. Finally, if no point of $`\mathrm{Spec}_{\text{ess}}(A)`$ lies strictly to the left of $`\lambda `$, then $`\lambda _{}=\mathrm{}`$ and there exist eigenvalues $`\lambda _j<\lambda `$ of $`A`$ converging to $`\mathrm{}`$. Here we set $`_j^{}:=E(\{\lambda _j\})`$. In all the above cases $`_j^{}E((\mathrm{},\lambda ))`$, $`j`$.

Similarly one can define subspaces $`_j^+E((\lambda ,+\mathrm{}))`$, $`j`$. It is clear that $`_j^\pm `$ are all mutually orthogonal, $`_j^\pm \mathrm{Dom}(A)`$, $`A_j^\pm _j^\pm `$, and there exist $`r_k>0`$ such that $`r_k+\mathrm{}`$ and $`_j^\pm E([r_k,r_k])`$, $`jk`$. Let $`_0`$ be the closed linear subspace spanned by $`_j^\pm `$, $`j`$ and $`_0^{}`$ be its orthogonal complement.

Let $`v_k^\pm _k^\pm `$, $`v_k^\pm =1`$ and $`\mu _k^\pm :=(Av_k^\pm ,v_k^\pm )`$. Then it follows from the spectral theorem that $`\mu _k^{}<\lambda <\mu _k^+`$, $`k`$. Let $`_k`$ be the closed linear subspace spanned by $`_j^\pm `$, $`j<k`$, $`E([r_k,r_k])_0^{}`$, and the vector

$$u_k:=\sqrt{\frac{\mu _k^+\lambda }{\mu _k^+\mu _k^{}}}v_k^{}+\sqrt{\frac{\lambda \mu _k^{}}{\mu _k^+\mu _k^{}}}v_k^+.$$

It is easily seen that the sequence $`(_k)_k`$ is increasing and belongs to $`\mathrm{\Lambda }(A)`$. Since $`v_k^\pm ,Av_k^\pm _k^\pm E([r_k,r_k])`$, the orthogonality properties of $`_j^\pm `$ imply

$`u_k=1,`$
$`((A\lambda I)u_k,u_k)={\displaystyle \frac{\mu _k^+\lambda }{\mu _k^+\mu _k^{}}}(Av_k^{},v_k^{})+{\displaystyle \frac{\lambda \mu _k^{}}{\mu _k^+\mu _k^{}}}(Av_k^+,v_k^+)`$
$`\lambda u_k^2={\displaystyle \frac{\mu _k^+\lambda }{\mu _k^+\mu _k^{}}}\mu _k^{}+{\displaystyle \frac{\lambda \mu _k^{}}{\mu _k^+\mu _k^{}}}\mu _k^+\lambda =0,`$
$`((A\lambda I)u_k,v)=0,vE([r_k,r_k])_0^{}{\displaystyle \underset{j<k}{}}(_j^{}_j^+).`$
Hence $`((A\lambda I)u_k,u)=0`$, $`u_k`$, i.e. $`P_k(A\lambda I)u_k=0`$, i.e. $`P_kAu_k=\lambda u_k`$, i.e. $`\lambda `$ is an eigenvalue of $`P_kA:_k_k`$, $`k`$. ∎

###### Remark 5.1.

If $``$ is separable, the subspaces $`_k`$ in Theorem 2.1 can be made finite dimensional. Indeed, let $`e_{j,m}^\pm `$, $`m=1,\mathrm{}`$ be an orthonormal basis of $`_j^\pm `$ and $`e_m^0`$, $`m=1,\mathrm{}`$ be an orthonormal basis of $`_0^{}`$. These bases exist because all subspaces of $``$ are separable. Let $`_k`$ be the closed linear subspace spanned by $`e_{j,m}^\pm `$ and $`E([r_k,r_k])e_m^0`$ with $`j,m<k`$, and the vector

$$u_k:=\sqrt{\frac{\mu _k^+\lambda }{\mu _k^+\mu _k^{}}}e_{k,1}^{}+\sqrt{\frac{\lambda \mu _k^{}}{\mu _k^+\mu _k^{}}}e_{k,1}^+.$$

Then it is clear that $`_k`$ is finite dimensional, and the rest of the proof of Theorem 2.1 remains unchanged.

###### Proof of Theorem 2.2.

It is sufficient to prove that for any finite dimensional $`\mathrm{Dom}(A)`$ and any $`\epsilon ,R>0`$ there exists $`^{}\mathrm{Dom}(A)`$ such that

$$PP^{}<\epsilon \text{and}\mathrm{Spec}(A,^{})(R,+\mathrm{}),$$

where $`P:`$ and $`P^{}:^{}`$ are the corresponding orthogonal projections.

Let $`m:=dim`$ and $`\phi _1,\mathrm{},\phi _m`$ be an orthonormal basis of $``$. Since $`A`$ is not bounded above, for any $`\tau >0`$ there exist mutually orthogonal vectors $`\psi _1,\mathrm{},\psi _m\mathrm{Dom}(A)E((\tau ,+\mathrm{}))`$ which we normalize by $`\psi _l=\delta `$, $`l=1,\mathrm{},m`$, where $`\delta (0,1)`$. The numbers $`\tau `$ and $`\delta `$ will be chosen later.

It is easy to see that the vectors $`\phi _1+\psi _1,\mathrm{},\phi _m+\psi _m`$ are linearly independent. Indeed, for any $`c_1,\mathrm{},c_m`$ Parseval’s identity implies

$$\underset{l=1}{\overset{m}{}}c_l(\phi _l+\psi _l)\underset{l=1}{\overset{m}{}}c_l\phi _l\underset{l=1}{\overset{m}{}}c_l\psi _l=(1\delta )\left(\underset{l=1}{\overset{m}{}}|c_l|^2\right)^{1/2}.$$
(10)
Similarly

$$\underset{l=1}{\overset{m}{}}c_l(\phi _l+\psi _l)(1+\delta )\left(\underset{l=1}{\overset{m}{}}|c_l|^2\right)^{1/2}.$$
(11)
Let $`^{}`$ be the subspace spanned by $`\phi _1+\psi _1,\mathrm{},\phi _m+\psi _m`$ and consider the operators $`J:^{}`$ and $`J^{}:^{}`$ defined by

$$J\left(\underset{l=1}{\overset{m}{}}c_l\phi _l\right)=\underset{l=1}{\overset{m}{}}c_l(\phi _l+\psi _l),J^{}\left(\underset{l=1}{\overset{m}{}}c_l(\phi _l+\psi _l)\right)=\underset{l=1}{\overset{m}{}}c_l\phi _l.$$

It is clear that $`J^{}J=I_{}`$, $`JJ^{}=I_{^{}}`$ and

$$(IJ)\left(\underset{l=1}{\overset{m}{}}c_l\phi _l\right)=\underset{l=1}{\overset{m}{}}c_l\psi _l,(IJ^{})\left(\underset{l=1}{\overset{m}{}}c_l(\phi _l+\psi _l)\right)=\underset{l=1}{\overset{m}{}}c_l\psi _l.$$

It follows from Parseval’s identity and (10), (11) that

$$J1+\delta ,J^{}\frac{1}{1\delta },(IJ)P\delta ,(IJ^{})P^{}\frac{\delta }{1\delta }.$$

Further,

$$PP^{}PP^{}P+P^{}PP^{}=PP^{}P+PP^{}P^{},$$

since $`(P^{}PP^{})^{}=PP^{}P^{}`$. Further, for any $`x`$ the closest point of $`^{}`$ is $`P^{}x`$. Hence $`PP^{}PPJP`$. Similarly $`P^{}PP^{}P^{}J^{}P^{}`$. Consequently

$$PP^{}(IJ)P+(IJ^{})P^{}\delta +\frac{\delta }{1\delta }=\frac{\delta (2\delta )}{1\delta }.$$

Fix $`\delta (0,1)`$ such that $`\delta (2\delta )/(1\delta )<\epsilon `$. Then $`PP^{}<\epsilon `$.

Since $``$ is a finite dimensional subspace of $`\mathrm{Dom}(A)`$, there exists $`M>0`$ such that $`AuMu`$, $`u`$. Take an arbitrary

$$v=\underset{l=1}{\overset{m}{}}c_l(\phi _l+\psi _l)^{}$$

and set

$$\phi :=\underset{l=1}{\overset{m}{}}c_l\phi _l,\psi :=\underset{l=1}{\overset{m}{}}c_l\psi _l\mathrm{Dom}(A)E((\tau ,+\mathrm{})).$$

Then

$`(Av,v)=(A\psi ,\psi )+2\mathrm{Re}(A\phi ,\psi )+(A\phi ,\phi )`$
$`\tau \psi ^22M\phi \psi M\phi ^2=(\tau \delta ^22M\delta M){\displaystyle \underset{l=1}{\overset{m}{}}}|c_l|^2.`$
Choose $`\tau `$ such that $`(\tau \delta ^22M\delta M)/(1+\delta )^2R+1`$. Then

$$(Av,v)\frac{\tau \delta ^22M\delta M}{(1+\delta )^2}v^2(R+1)v^2$$

(see (11)). Hence $`\mathrm{Spec}(A,^{})[R+1,+\mathrm{})`$. ∎

Suppose $`A`$ is a bounded below self-adjoint operator. Then everything said above generalizes from $`\mathrm{\Lambda }(A)`$ to $`\mathrm{\Lambda }(Q_A)`$ (see Section 2 for the definition) with almost no changes. One only needs to modify the concluding estimates of the proof of Theorem 2.2 in the following way. Adding if necessary a sufficiently large positive number to the bounded below operator $`A`$, one can assume without loss of generality that $`A`$ is positive definite. Then

$`Q_A(v,v)=Q_A(\psi ,\psi )+2\mathrm{Re}Q_A(\phi ,\psi )+Q_A(\phi ,\phi )`$
$`Q_A(\psi ,\psi )2Q_A(\psi ,\psi )^{1/2}Q_A(\phi ,\phi )^{1/2}`$
$`{\displaystyle \frac{1}{2}}Q_A(\psi ,\psi )2Q_A(\phi ,\phi ){\displaystyle \frac{\tau }{2}}\psi ^22M\phi ^2=`$
$`{\displaystyle \frac{\tau \delta ^24M}{2}}{\displaystyle \underset{l=1}{\overset{m}{}}}|c_l|^2.`$
Choose $`\tau `$ such that $`(\tau \delta ^24M)/2(1+\delta )^2R+1`$. Then

$$Q_A(v,v)\frac{\tau \delta ^24M}{2(1+\delta )^2}v^2(R+1)v^2.$$

Hence $`\mathrm{Spec}(A,^{})[R+1,+\mathrm{})`$.

Now we proceed to the proof of Theorem 2.5. For any $`c_1,c_2`$ we denote by $`𝔻(c_1,c_2)`$ the open disk with diameter $`(c_1,c_2)`$.

###### Lemma 5.2.

Let $`A`$ be a self-adjoint operator on a Hilbert space $``$ and let $``$ be a finite dimensional subspace of $`\mathrm{Dom}(A)`$. Suppose $`a,b`$ are such that

$$(a,b)\mathrm{Spec}(A)=\mathrm{}.$$

Then

$$𝔻(a,b)\mathrm{Spec}_2(A,)=\mathrm{}.$$

###### Proof of Lemma 5.2.

Take an arbitrary $`z𝔻(a,b)`$. The interval $`(a,b)`$ is seen from $`z`$ at an angle $`\pi /2+\epsilon `$, where $`\epsilon (0,\pi /2]`$. It is clear that the set

$$F=\{(\lambda z)^2:\lambda \mathrm{Spec}(A)\}$$

lies in a sector with vertex at 0 and angle $`\pi 2\epsilon `$. The distance from 0 to $`F`$ is obviously positive. Therefore there exist $`\vartheta (\pi ,\pi ]`$ and $`c>0`$ such that

$$\mathrm{Re}\left(e^{i\vartheta }(\lambda z)^2\right)c,\lambda \mathrm{Spec}(A).$$

Now applying the spectral theorem for $`A`$ one obtains

$`\mathrm{Re}\left(e^{i\vartheta }((AzI)u,(A\overline{z}I)u)\right)=`$
$`\mathrm{Re}\left(e^{i\vartheta }\left(Au^22z(Au,u)+z^2u^2\right)\right)=`$
$`\mathrm{Re}\left(e^{i\vartheta }{\displaystyle _{\text{Spec}(A)}}(\lambda ^22\lambda z+z^2)d(E(\lambda )u,u)\right)=`$
$`{\displaystyle _{\text{Spec}(A)}}\mathrm{Re}\left(e^{i\vartheta }(\lambda z)^2\right)d(E(\lambda )u,u)c{\displaystyle _{\text{Spec}(A)}}d(E(\lambda )u,u)=`$
$`cu^2,u.`$
Hence $`u\{0\}`$ cannot satisfy (4), i.e. $`z\mathrm{Spec}_2(A,)`$. ∎

###### Proof of Theorem 2.5.

Suppose the intersection in the statement of the Theorem is empty. Then there exists $`\epsilon >0`$ such that

$$\mathrm{Spec}(A)[\mathrm{Re}z|\mathrm{Im}z|\epsilon ,\mathrm{Re}z+|\mathrm{Im}z|+\epsilon ]=\mathrm{}.$$

Hence $`𝔻(\mathrm{Re}z|\mathrm{Im}z|\epsilon ,\mathrm{Re}z+|\mathrm{Im}z|+\epsilon )`$ does not intersect $`\mathrm{Spec}_2(A,)`$ according to Lemma 5.2. On the other hand, $`z`$ belongs to both of these sets. The obtained contradiction proves the theorem. ∎

## 6 More on Example I

Here we describe what happens if the function $`a`$ of Section 3 is slightly modified. Since none of the changes we consider below affects the second order relative spectra in any serious way, we concentrate on the standard (first order) projection method.

We have performed calculations, similar to those presented in Section 3, for the operator $`A`$ corresponding to the function which is obtained from $`a`$ by multiplying the cosines by 0.01, 0.001, or zero. The only noticeable difference with the results in Section 3 is that this time there are more spurious eigenvalues in the gap of the essential spectrum and additional spurious eigenvalues appear, as $`N`$ increases, quicker than in the original example (cf. Table 1 and Remark 3.1).

The situation becomes more interesting if one moves the discontinuity of $`a`$ from 0 to a point $`y0`$:

$$a(x)=a_y(x)=\{\begin{array}{cc}\frac{3}{2}+\frac{1}{2}\mathrm{cos}\sqrt{5}x,\hfill & \text{for}\pi x<y,\hfill \\ 2+\mathrm{cos}\sqrt{2}x,\hfill & \text{for}yx<\pi .\hfill \end{array}$$
(12)
The corresponding operator of multiplication by $`a_y`$ will be denoted $`A_y`$.

To put the results into the right perspective, we need to return to the original example $`A=A_0`$ first. We know from Section 3 that $`\mathrm{Spec}(A,_N)`$ show remarkable stability as $`N=2n+1`$ increases. The first surprise happens when one looks at $`\mathrm{Spec}(A,_N)`$ for $`N=2n`$, where $`_N:=\text{span}(\{e_k\}_{k=n}^{n1})`$. These also seem to be converging, but the spurious eigenvalues differ considerably from those corresponding to $`N=2n+1`$ (see Figure 12 and Table 2 and compare them with Figure 2 and Table 1).

One might argue that this difference between the two cases occurs due to the nonsymmetric truncation in the case $`N=2n`$ and that taking only odd values of $`N`$ is more natural. We prefer, however, to interpret this fact as $`2`$-periodicity in $`N`$ of the points of $`\mathrm{Spec}(A,_N)`$ lying in the gap $`(1,1)`$ of the essential spectrum of $`A`$, i.e., of the spurious eigenvalues of $`A`$. Such an interpretation is supported by Table 3, which shows the spurious eigenvalues for some consecutive values of $`N`$.

We show below that in a more general case of the operator $`A_y`$ the “period” depends on the location of the point $`y`$ of the discontinuity of $`a_y`$ (see (12)). This period is equal to $`2`$ in the original example only because we take $`y=0`$ there.

Let $`y=\pi /3`$. Figure 13 shows the computed (spurious) eigenvalues of $`A_y=a_yI`$ in the interval $`(1,1)`$ as functions of $`N`$. Remarkably, the eigenvalues vary quite slowly with $`N`$ as long as $`N=\text{const}(mod3)`$. To emphasize this “3–periodicity” in $`N`$, we use different symbols for points with different values of $`N(mod3)`$.

A similar effect occurs for the operator $`B_y`$, which is obtained from the operator $`B`$ in (5) by replacing $`a`$ with $`a_y`$ (see Figure 14). However, curves corresponding to the eigenvalues in the interval $`(1,1)`$, which were confirmed by the second order relative spectra as the correct ones, do not depend on the value of $`N(mod3)`$. They also seem to converge to constants faster than those depending on $`N(mod3)`$.

We have performed similar calculations with $`y=\pi /(p/q)=\pi q/p`$ for the following values $`p/q=`$ 2, 4, 5, 6, 7, 21, 21/10, 21/5, 7/3, 8/3, and 5/2 (see Figure 15 for the results for some of the values). In all these cases the spurious eigenvalues of $`A_y`$ showed a very clear pattern of “periodic” behaviour with the minimal “period”

$$\omega (q/p):=\{\begin{array}{cc}2,\hfill & \text{if}q=0,\hfill \\ p,\hfill & \text{if}p\text{and}q\text{are both odd,}\hfill \\ 2p,\hfill & \text{if}q0\text{and either}p\text{or}q\text{is even.}\hfill \end{array}$$

Such “periodic” behaviour can be observed in Figure 15. The appearance of $`\omega (q/p)`$ here is probably not very surprising, as it is, for any simple fraction $`q/p`$, the smallest natural number such that $`\omega \times `$(length of $`[\pi ,\pi q/p]`$) and $`\omega \times `$(length of $`[\pi q/p,\pi ]`$) are integer multiples of $`2\pi `$.

## 7 More on Example II

As we have shown in Section 4, a FEM calculation of the spectrum of operator $`H`$ using the projection space $`𝐋_N\times 𝐋_N`$ leads to heavy spectral pollution in the interval $`(1,2)`$. At the same time, the second order spectra relative to the projection space $`𝐌_N\times 𝐋_N`$ approximate the spectrum of $`H`$ very well.

A natural question is whether such a discrepancy may be mostly due to the different projection spaces being used in the calculations of the first order and second order spectra. To check this, we compute the spectra $`\mathrm{Spec}(H,𝐌_N\times 𝐋_N)`$ using the same projection spaces as we used in Section 4 for the second order relative spectra. The results (restricted to eigenvalues lying in the interval $`[0,6]`$) are shown in Figure 16. This time, there are no spurious eigenvalues at all! In the particular case $`N=100`$, we have 100 eigenvalues in $`[0,1]`$, and 200 eigenvalues in $`[2,\mathrm{})`$ (recall that $`dim(𝐌_N\times 𝐋_N)=3N`$). The actual eigenvalues of $`H`$ are well approximated in this case as well, see Remark 7.1 below.

This of course indicates that the projection spaces $`𝐌_N\times 𝐋_N`$ are much better suited for the calculation of the spectrum of $`H`$ than the projection spaces $`𝐋_N\times 𝐋_N`$. Such a dependence upon the choice of the spaces of finite elements is quite interesting, especially because it does not occur for the operator $`d^2/dx^2`$ and because of the explicit relations which exist between the eigenvalues and eigenfunctions of $`H`$ and those of $`d^2/dx^2`$, see (8) and (9).

The spaces $`𝐌_N\times 𝐋_N`$ are well suited for the calculation of $`\mathrm{Spec}(H)`$ because they are subspaces of the mixed order Sobolev space $`W_2^2([\pi ,\pi )])\times W_2^1([\pi ,\pi )])`$ and therefore reflect much better the mixed orders of differentiation in the operator $`H`$, see also \[RSHSPV\] and references therein. We attempt to explain this phenomenon in a more general setting below, see Remark 7.2.

One may be tempted to think that the success of the second order method used in Section 4 is also only due to the correct choice of the finite elements. In order to check this assumption let us compare the results of calculation of the first order and second order relative spectra using the projection spaces $`𝐌_N\times 𝐌_N`$. The results for $`\mathrm{Spec}(H,𝐌_N\times 𝐌_N)`$ are shown in Figure 17 — again, there is spectral pollution in the interval $`(1,2)`$.

On the other hand, the calculations of the second order relative spectra $`\mathrm{Spec}_2(H,𝐌_N\times 𝐌_N)`$ produce pictures very similar to those in Figures 10 and 11. The only noticeable difference is the presence of a number of points with relatively large moduli of imaginary parts ($`10`$ for $`N=50`$) which seem to lie on a “nice” curve passing through the point $`\lambda =1`$ (see Figure 18; note the different scale along the imaginary axis compared with either Figure 10 or 11).

Thus, for the spaces $`𝐌_N\times 𝐌_N`$, the projection method using second order relative spectra works quite well, whereas the standard method again fails dramatically.

###### Remark 7.1.

It should be noted that both the first order and the second order relative spectra (outside $`(1,2)`$) approximate $`\mathrm{Spec}(H)`$ very well, with accuracy depending on the finite elements chosen. The points of $`\mathrm{Spec}(H,𝐋_N\times 𝐋_N)`$, $`\mathrm{Spec}(H,𝐌_N\times 𝐋_N)`$, and $`\mathrm{Spec}(H,𝐌_N\times 𝐌_N)`$ converge to the points $`\lambda _k^\pm \mathrm{Spec}(H)`$ with the rate of about $`O(N^2)`$, $`O(N^4)`$, and $`O(N^6)`$, respectively. For the second and third choice of finite elements, the same rates of convergence are achieved by the real parts of the points of $`\mathrm{Spec}_2(H,𝐌_N\times 𝐋_N)`$ and $`\mathrm{Spec}_2(H,𝐌_N\times 𝐌_N)`$, whereas the imaginary parts converge to zero with the rate of about $`O(N^2)`$ and $`O(N^3)`$, respectively. These convergence rates were obtained by analyzing the $`\mathrm{log}`$$`\mathrm{log}`$ dependence of errors upon $`N`$ for several low eigenvalues taken from both $`\lambda ^+`$ and $`\lambda ^{}`$ series.

###### Remark 7.2.

Here we try to explain why the standard projection method of approximating $`\mathrm{Spec}(H)`$ does not lead to spectral pollution in the case of the projection spaces $`𝐌_N\times 𝐋_N`$. The main reason is that it seems, in this particular case, to “mimic” the way one finds $`\mathrm{Spec}(H)`$ analytically. Indeed, let $`_1`$ be a space of finite elements for $`v`$, $`_2`$ be a space of finite elements for $`u`$, and let $`_2^{}:=\{\psi ^{}:\psi _2\}`$. Let $`\phi _j`$, $`j=1,\mathrm{},M`$ be a basis of $`_1`$ and $`\psi _k`$, $`k=1,\mathrm{},N`$ be a basis of $`_2`$. Since the projection method does not depend on the choice of bases, provided all the calculations are precise, we can assume that the basis $`\{\phi _j\}`$ is orthonormal. Then the finite dimensional system for (6) takes the form

$$\{\begin{array}{cc}\hfill _{k=1}^N(\psi _m^{},\psi _k^{})u_k+_{j=1}^M(\psi _m^{},\phi _j)v_j=& \lambda _{k=1}^N(\psi _m,\psi _k)u_k,\hfill \\ & m=1,\mathrm{},N,\hfill \\ \hfill _{k=1}^N(\phi _j,\psi _k^{})u_k+2v_j=& \lambda v_j,\hfill \\ & j=1,\mathrm{},M.\hfill \end{array}$$

This leads to the following system for the coefficients $`u_k`$:

$`{\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m^{},\psi _k^{})u_k+{\displaystyle \frac{1}{\lambda 2}}{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{M}{}}}(\psi _m^{},\phi _j)(\phi _j,\psi _k^{})u_k=`$
$`\lambda {\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m,\psi _k)u_k,m=1,\mathrm{},N.`$ (13)
Let $`P_1`$ be the orthogonal projection onto $`_1`$. Then

$$P_1\psi _k^{}=\underset{j=1}{\overset{M}{}}(\psi _k^{},\phi _j)\phi _j\text{and}(P_1\psi _m^{},P_1\psi _k^{})=\underset{j=1}{\overset{M}{}}(\psi _m^{},\phi _j)(\phi _j,\psi _k^{}).$$

Hence (7.2) is equivalent to

$`{\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m^{},\psi _k^{})u_k+{\displaystyle \frac{1}{\lambda 2}}{\displaystyle \underset{k=1}{\overset{N}{}}}(P_1\psi _m^{},P_1\psi _k^{})u_k=\lambda {\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m,\psi _k)u_k,`$ (14)
$`m=1,\mathrm{},N.`$
On the other hand, it is clear that (6) implies

$$u^{\prime \prime }\frac{1}{\lambda 2}u^{\prime \prime }=\lambda u,$$
(15)
which is equivalent to (7). So, the projection method applied to (15) should not lead to spectral pollution. This gives the following finite dimensional system

$`{\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m^{},\psi _k^{})u_k+{\displaystyle \frac{1}{\lambda 2}}{\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m^{},\psi _k^{})u_k=\lambda {\displaystyle \underset{k=1}{\overset{N}{}}}(\psi _m,\psi _k)u_k,`$ (16)
$`m=1,\mathrm{},N.`$
If $`_1_2^{}`$, then (14) and (16) coincide and the projection method applied to (6) does not pollute. If $`_1_2^{}`$, but $`_1`$ and $`_2^{}`$ are “compatible”, then it may well happen that (14) mimics the properties of (16) at least partially, and the projection method applied to (6) does not lead to spectral pollution.

## 8 Concluding remarks

The main message of the paper is that projection methods using second order relative spectra may provide an efficient tool for dealing with spectral pollution. They are not much harder to implement than the standard projection methods and are likely to give useful a posteriori information (see Theorem 2.5) about the location of (at least a part of) the spectrum of a self-adjoint operator. We would feel that the paper has achieved its goal if it persuaded other people to try second order relative spectra in their problems of interest. We hope that the method we suggest here will work well in problems more serious than the model examples considered above.

We would like to stress again that we do not claim that second order relative spectra are the only way of addressing spectral pollution. In our Example I, for instance, one could argue that spectral pollution could have been detected just by looking at pictures like Figure 3, without any reference to the second order relative spectra. The argument would be that a genuine eigenfunction “cannot” look like the one in Figure 3, and hence the corresponding eigenvalue is spurious. The “periodicity” of spurious eigenvalues (see Section 6) might also have been used to detect possible spectral pollution. Indeed, those approximate eigenvalues of $`B_y`$, which were confirmed by the second order relative spectra as the correct ones, do not exhibit “periodicity” in $`N`$. They are present for any sufficiently large value of $`N`$. We have no intention to disagree with these arguments. Note however that both of these approaches are based on the results of calculations using the standard (first order) projection method, which does not see the difference between $`A`$ and the Toeplitz operator $`T_a`$. Consequently, they cannot distinguish between spectral pollution for the former and non-pollution for the latter. Second order relative spectra, on the other hand, are free from this problem. Our point is that it is very easy to fail to notice spectral pollution unless one is very careful, and that second order relative spectra seem to be very efficient “lie detectors”.

Another perfectly legitimate argument is that one should try to avoid spectral pollution by making a sensible choice of the subspaces $`_k`$, rather than allow it to happen and then try to detect it. The standard projection method using $`𝐌_N\times 𝐋_N`$, for instance, works very well in our Example II (see Section 7). One should, however, bear in mind that choosing “right” subspaces when applying the standard projection method in more difficult situations may prove non-trivial and require a good understanding of the problem at hand. Second order relative spectra, on the other hand, seem to provide a more robust projection method which works well even in the case of “wrong” subspaces, when the standard one fails. The case of $`𝐌_N\times 𝐌_N`$ considered in Section 7 is an example of such situation.

Finally, we refer the reader to the companion paper by E. B. Davies and M. Plum \[DP\], where a different method of treating spectral pollution is developed.

## Acknowledgements

We would like to thank Brian Davies for interest in our work and very helpful discussions. This work was partially supported by a grant from the EPSRC Spectral Theory Network.