High-statistics study of K^0_S pair production in two-photon collisions

We report a high-statistics measurement of the differential cross section of the process gamma gamma -->K^0_S K^0_S in the range 1.05 GeV<= W<= 4.00 GeV, where W is the center-of-mass energy of the colliding photons, using 972 fb^{-1} of data collected with the Belle detector at the KEKB asymmetric-energy e^+ e^- collider operated at and near the Upsilon-resonance region. The differential cross section is fitted by parameterized S-, D_0-, D_2-, G_0- and G_2-wave amplitudes. In the D_2 wave, the f_2(1270), a_2(1320) and f_2'(1525) are dominant and a resonance, the f_2(2200), is also present. The f_0(1710) and possibly the f_0(2500) are seen in the S wave. The mass, total width and product of the two-photon partial decay width and decay branching fraction to the K bar{K} state Gamma_{gamma gamma}B(K bar{K}) are extracted for the f_2'(1525), f_0(1710), f_2(2200) and f_0(2500). The destructive interference between the f_2(1270) and a_2(1320) is confirmed by measuring their relative phase. The parameters of the charmonium states chi_{c0} and chi_{c2} are updated. Possible contributions from the chi_{c0}(2P) and chi_{c2}(2P) states are discussed. A new upper limit for the branching fraction of the P- and CP-violating decay channel eta_c -->K^0_S K^0_S is reported. The detailed behavior of the cross section is updated and compared with QCD-based calculations.

We report a high-statistics measurement of the differential cross section of the process γγ → K 0 S K 0 S in the range 1.05 GeV ≤ W ≤ 4.00 GeV, where W is the center-of-mass energy of the colliding photons, using 972 fb −1 of data collected with the Belle detector at the KEKB asymmetric-energy e + e − collider operated at and near the Υ-resonance region. The differential cross section is fitted by parameterized S-, D0-, D2-, G0-and G2-wave amplitudes. In the D2 wave, the f2(1270), a2(1320) and f 2 (1525) are dominant and a resonance, the f2(2200), is also present. The f0(1710) and possibly the f0(2500) are seen in the S wave. The mass, total width and product of the two-photon partial decay width and decay branching fraction to the KK state ΓγγB(KK) are extracted for the f 2 (1525), f0(1710), f2(2200) and f0(2500). The destructive interference between the f2(1270) and a2(1320) is confirmed by measuring their relative phase. The parameters of the charmonium states χc0 and χc2 are updated. Possible contributions from the χc0(2P ) and χc2(2P ) states are discussed. A new upper limit for the branching fraction of the P -and CP -violating decay channel ηc → K 0 S K 0 S is reported. The detailed behavior of the cross section is updated and compared with QCD-based calculations.

I. INTRODUCTION
We present a high-statistics study of the cross section for the process γγ → K 0 S K 0 S , through the measurement of e + e − → (e + e − )K 0 S K 0 S where neither a scattered electron nor positron is detected (zero-tag mode), in the W region from close to its threshold to 4.0 GeV and in the angular range | cos θ * | ≤ 0.8, where W is the total energy of the parent photons and θ * is the scattering angle of the K 0 S in their center-of-mass (c.m.) reference frame. Measurements of exclusive hadronic final states in twophoton collisions provide valuable information concerning the physics of light-and heavy-quark resonances, perturbative and non-perturbative QCD and hadron-production mechanisms. The Belle collaboration has measured the production cross sections for charged-pion pairs [1][2][3], charged and neutral-kaon pairs [3][4][5], and proton-antiproton pairs [6]. Belle has also analyzed D-meson-pair production and observed a new charmonium state identified as the χ c2 (2P ) [7]. In addition, Belle has measured the production cross section for the π 0 π 0 , ηπ 0 and ηη final states [8][9][10][11]. The statistics of these measurements are two to three orders of magnitude higher than in pre-B-factory measurements [12], opening a new era in two-photon physics.
The f J and a J mesons (with even spin J) both contribute to the process of γγ → KK. The almost degenerate f J and a J that are predominantly uū and dd are predicted to interfere destructively in γγ → K 0K 0 and constructively in γγ → K + K − [13]. This is due to the Okubo-Zweig-Iizuka rule [14] where the dd (uū) initial state dominates in K 0K 0 (K + K − ) production. To the extent that the ss component is ignored, the dd (uū) state can be expressed as (f J − a J )/ √ 2 ((f J + a J )/ √ 2) by the isospin consideration. In the γγ → K 0 S K 0 S reaction near the threshold, Refs. [15,16] predict a destructive interference between the f 0 (980) and a 0 (980), irrespective of their nature, that suppresses the production cross section to below 1 nb. They consider the K 0 S K 0 S production to be dominated by the rescattering process of K + K − → K 0K 0 near the threshold. There have been no further data to shed light on this.
The third was attributed to the f J (1710) (J = 2) [21]. The limited statistics of these experiments (e.g., 0.588 fb −1 for the L3 results [21]) were insufficient to resolve and to study higher mass resonances. Although these experiments operated at higher e + e − c.m. energies, the cross section of each two-photon production process in a specific W range rises only logarithmically with the e + e − c.m. energy. The CLEO collaboration published the distribution of the invariant mass for γγ → K 0 S K 0 S in a search for η(1440) → K 0 S K ± π ∓ based on 13.8 fb −1 of data [22]; the K 0 S K 0 S measurement was used solely for the calibration of the K 0 S efficiency, but no physics results were extracted. Intriguingly, several resonant structures can be observed clearly in their K 0 S K 0 S mass spectrum. In the previous Belle study of the γγ → K + K − reaction, enhancements near 1.75 GeV, 2.0 GeV and 2.3 GeV were reported and attributed to the a 2 (1700), f 2 (2010) and f 2 (2300), respectively [4,23].
In this article, we present a high-statistics study of the cross section for γγ → K 0 S K 0 S from close to its threshold to W = 4.0 GeV. The data are based on an integrated luminosity of 972 fb −1 . This significantly extends our previous study [5], where the measurement of this process was reported for 2.4 GeV ≤ W ≤ 4.0 GeV with an integrated luminosity of 397.6 fb −1 . In that study, we compared the high-energy behavior of the cross section with the QCD-based calculations or models [24,25]. Signals for the χ c0 and χ c2 charmonium states were observed. Here, we extend the c.m. energy lower limit down to 1.05 GeV and investigate the intermediatemass resonances with higher statistics data.
We report the first measurement of the differential cross section for γγ → K 0 S K 0 S below 2.4 GeV. Previously, only the event distributions were obtained for this process [17][18][19]21] and the integrated cross section was presented with limited statistics [20]. In analyzing the differential cross section, we measure the phase difference between the a 2 (1320) and f 2 (1270) as well as the parameters (mass, width and product of the two-photon partial decay width and decay branching fraction to the KK, Γ γγ B(KK) ) of the f 2 (1525) including the interference. Resonance-like enhancements are investigated in the region W > 1.6 GeV. We also provide some new information on possible glueball candidates such as the f 0 (1710) and f J (2220) [26][27][28][29].
We then update the measurements of the parameters of the χ c0 and χ c2 states. Possible contributions from the radially excited states χ cJ (2P ) are investigated. The χ c2 (2P ) was discovered and confirmed in two-photon collisions [23], and the X(3915) found in the γγ → ωJ/ψ process has been identified re-cently as the χ c0 (2P ) state [30]. In addition, we also report searches for the P -and CP -violating decay η c → K 0 S K 0 S and set a new upper limit for its branching fraction. Finally, we compare the cross section dependence on W and | cos θ * | for W > 2. 6 GeV with QCD predictions.
This article is organized as follows. First we describe the details of the data selection (Sec. II), background subtraction (Sec. III), efficiency determination (Sec. IV) and derivation of the differential cross section (Sec. V). We then present results on resonance analysis (Sec. VI), update the properties of several charmonia (Sec. VII), and model the crosssection behavior for W > 2.6 GeV (Sec. VIII). Finally, we present a summary and draw conclusions (Sec. IX).

II. THE EXPERIMENTAL APPARATUS AND SELECTION OF SIGNAL CANDIDATES
In this section, we describe the Belle detector, data sample, triggers, Monte Carlo simulation program and selection of signal candidates.

A. Experimental apparatus
Data were collected with the Belle detector operated at the KEKB asymmetric-energy e + e − collider [31,32]. A comprehensive description of the Belle detector is given elsewhere [33,34]. In this paper we briefly discuss only those detector components that are essential for the described measurement. Charged tracks are reconstructed from hit information in the silicon vertex detector and the central drift chamber (CDC). The CDC is used as the main device to trigger readout for the events with charged particles. A barrel-like arrangement of time-of-flight (TOF) counters and trigger scintillation counters (TSC) are used to supplement the CDC trigger on charged particles and to measure their time of flight. Particle identification (ID) is achieved by including information from an array of aerogel threshold Cherenkov counters. Photon detection and energy measurements are performed with a CsI(Tl) electromagnetic calorimeter (ECL). All of the above detectors are located inside a superconducting solenoid coil that provides a uniform 1.5 T magnetic field. The detector solenoid is oriented along the z axis, pointing in the direction opposite that of the positron beam. The rϕ plane is transverse to this axis.

B. Data sample
This analysis is based on a data sample corresponding to an integrated e + e − luminosity of 972 fb −1 . Data were collected at the energy of the Υ(4S) resonance ( √ s = 10.58 GeV) and 60 MeV below it (784 fb −1 ), at energies between 10.6 GeV and 11.1 GeV (151 fb −1 , mainly near the Υ(5S) resonance at 10.88 GeV), and at lower energies between 9.4 GeV and 10.3 GeV (38 fb −1 , primarily near the Υ(2S) resonance at 10.02 GeV). We analyze these data with a common algorithm for selecting K 0 S pair candidates from a zero-tag two-photon process because the process is independent of incident e ± energies.

C. Triggers and filtering
The analysis is based on data recorded with triggers that are sensitive to low-transverse-momentum (p t ) pions from K 0 S → π + π − decays. Signal low-p t pions have large curvatures in the CDC and deposit only a small amount of energy in the ECL; as a result, the trigger efficiency for the signal pions decreases steeply toward the threshold energy for K 0 S K 0 S production. To reduce the uncertainty in the trigger efficiency, we select data events recorded inclusively with triggers A, B and C as described below. These triggers make use of full-(short-) length charged tracks in the CDC volume that have p t > 0.3 GeV/c (0.2 GeV/c < p t < 0.3 GeV/c) [35]. Trigger A requires two or more full-length tracks in the CDC wire layers with an opening angle of roughly 135 • or larger in the rϕ plane [36], and at least two TOF/TSCmodule hits [37] and energy deposit with more than 0.11 GeV in at least one ECL trigger segment. Trigger B requires two CDC tracks, of which at least one track is a full-length one, with the opening angle requirement of trigger A, as well as a low-energy threshold condition (LowE [38]) of 0.5 GeV for the ECL total energy. By design, there is a large redundancy between triggers A and B. Trigger C is a three-track trigger with TOF/TSC-module and ECL segment/energy requirements. This trigger is sensitive to short and full tracks, but must have hits in the TOF and ECL. Details of the efficiencies and correlations of the three triggers are discussed in Sec. IV B.
To be recorded, a candidate event must pass the level-4 software trigger (L4, see Ref. [39]), in which a fast track-finding program reconstructs one or more tracks with transverse momentum p t > 0.3 GeV/c, each satisfying the requirements on the point of closest approach of the track to the z axis of dr < 1 cm and |dz| < 4 cm, where dr and dz are the distances between this point and the interaction point (IP) in the rϕ plane and along the z direction, respectively.

D. Monte Carlo simulation
The signal Monte Carlo (MC) events for e + e − → e + e − K 0 S K 0 S are generated using the MC code TREPS [40] at 81 fixed W points between 1.0 and 4.1 GeV and isotropically in | cos θ * |. Variables with (without) the asterisk represent observables in the c.m. (laboratory) reference frame. As we cannot measure the γγ collision axis directly, in the measurement we approximate it by the e + e − -collision axis in the e + e − c.m. frame.
In our simulation, we use the experimental setup and background files for runs at √ s = 10.58 GeV. To study the dependence of the analysis on run conditions and beam energy, we have generated additional signal MC events at 14 W points for each of the different run periods at √ s = 10.58 GeV, and at 12 and 6 W points with √ s = 10.88 GeV and 10.02 GeV, respectively. We embed background hit patterns from random trigger data into MC events, thus taking into account the efficiency dependence on run conditions.
In the signal MC generator, the Q 2 max parameter, a maximum virtuality of the incident space-like photons is set to 1.0 GeV 2 . The form factor σ γγ (0, Q 2 ) = σ γγ (0, 0)/(1+Q 2 /W 2 ) 2 is assumed. This assumption does not affect the results of our analysis, because we select events with Q 2 ≈ | p * t | < 0.1 GeV/c, thus requiring Q 2 /W 2 to be much smaller than 1, where | p * t | is the transverse momentum of the γγ system in the e + e − c.m. frame. Although the maximum Q 2 value determined from the requirement of the nondetection range of the scattered electron/positron is about 2 GeV 2 , the | p * t | condition applied to data limits Q 2 more tightly to be less than ∼ 0.01 GeV 2 . The Q 2 max = 1 GeV 2 used in the MC is larger than this experimental limit, and in this case the choice of Q 2 max in the MC does not affect the final γγ-based cross section results; i.e., the Q 2 max value is included in the definitions of the luminosity function calculated by TREPS, as well as in the efficiency. As a result, their effects are compensated in the cross section derivation (see Eq. (6)).
A sample of 400,000 events is generated at each W point per experimental setup. These events are then processed through the detector and trigger simulations and reconstructed using the same algorithms as for the real data. The decay of the K 0 S meson is managed in the GEANT-based detector simulation [42].

E. Selection criteria
We select K 0 S K 0 S two-photon event candidates in which each K 0 S decays to π + π − and neither scattered lepton is detected, i.e., in the zero-tag mode. Such candidates are required to contain exactly four charged tracks with small total transverse momentum in which two pairs of oppositely charged tracks form K 0 S candidates with vertices significantly away from the IP.
In order to reduce the background contribution from e + e − annihilation processes, the sum of the absolute momenta of the four tracks must be less than 6 GeV/c and the total energy of all ECL clusters must be less than 6 GeV.
To reduce the systematic uncertainty arising from reconstruction efficiency, we use only good-quality tracks that have p t > 0.1 GeV/c, dr < 5 cm and |dz| < 5 cm. The vector sum of the transverse momenta of the four tracks | p t | must be less than 0.2 GeV/c, using the azimuthal direction of the tracks at their closest approach to the nominal IP on its curved trajectory in the magnetic field. Each of the four tracks has to be identified as a pion from the particle-ID detectors with a K/π likelihood ratio: L(K)/(L(K) + L(π)) < 0.8. The pion identification efficiency is larger than 99% for p < 0.6 GeV/c and 95% for p = 0.8 GeV/c. To further reduce the annihilation contribution, the invariant mass of the four tracks with the pion mass assignment is required to be less than 5 GeV/c 2 . To eliminate backgrounds that include π 0 mesons, we require that there be no π 0 candidates with p t > 0.1 GeV/c and χ 2 < 4 in the mass-constrained fit of the available two-photon combinations.
Each pair of tracks forming a K 0 S candidate must have a difference in z coordinates at their point of closest approach in the rϕ plane, |∆z|, satisfying |∆z| < (p K +1.6) cm, where p K is the K 0 S momentum in GeV/c. The momentum dependence here incorporates the effect of resolution in the vertex determination. The reconstructed invariant mass of the two pions, M ππ , should satisfy |M ππ − m K | < 20 MeV/c 2 , where m K is the nominal K 0 S mass. We require a unique assignment of the four pions as the decay products from the two K 0 S by rejecting events that have ambiguous combinations. We further require that exactly two K 0 S candidates that are reconstructed from non-overlapping combinations of two charged tracks are found in the event. Figure 1 shows a two-dimensional plot of the two measured K 0 S masses where K1 and K2 are randomly assigned in each event.
To further reduce the background contribution 6 and to select well-reconstructed events, we require the difference of the reconstructed masses of the two K 0 S to satisfy |M K1 − M K2 | < 10 MeV/c 2 . We define the average of the reconstructed masses of the two These selection criteria are depicted in Fig. 1 with diagonal lines. Then, the decay position and momentum vector of each K 0 S are determined by a kinematical fit.
The radial displacement of each K 0 S vertex from the nominal IP, r V , must satisfy the condition r V i > max(0, W − 2) × 0.1 cm, where W is in GeV. This requirement does not apply to events with W < 2 GeV.
Backgrounds from the non-K 0 S K 0 S two-photon four-charged-pion production process (the "fourpion" process) are strongly suppressed if we require the two K 0 S vertices to be spatially separated, using combinations of two-dimensional (d V r ) and threedimensional (d V ) distances. The signed distance between the two vertices in the rϕ plane, d V r , defined according to must satisfy d V r > +0.05 cm, where r V i and p ti are two-dimensional vectors projected onto the rϕ plane of the decay vertex and transverse momentum, respectively, for each K 0 S . The event must satisfy either d V > 0.7 cm or d V r > +0.3 cm, where d V is a distance between the two vertices in the threedimensional space. Figures 2 and 3 show the distributions for these distances in the data (before the above selection criteria on them are applied) and signal MC samples. The peaks near zero in the data are due to the fourpion process γγ → π + π − π + π − whose cross section is larger than the signal one. This process is discussed in Secs. III B and IV B 3. Note that events with d V r < +0.05 cm or d V < 0.3 cm are rejected by our selection criteria and the relation |d V r | ≤ d V . We further require the projection of the distance between the vertices in the rϕ plane onto the vector of the transverse momentum difference, δ V , defined by to satisfy δ V < 0.7 cm, where ∆ϕ is the azimuthalangle difference between the vertex-position difference vector and the transverse-momentum difference vector.
To further eliminate events with significant photon activity, we require the total energy deposit in the  ECL to satisfy where E Ki is the total energy of each K 0 S . This selection criterion is determined by a study based on the signal MC in order not to lose any significant efficiency even if a pion deposits energy in the ECL after a nuclear interaction.
Finally, the p t balance of the K 0 S pair in the e + e − c.m. frame is required to satisfy | p * t | < 0.1 GeV/c.
We select candidates in the region 1.05 GeV ≤ W ≤ 4.10 GeV and | cos θ * | < 0.8. The W distribution of the selected K 0 S K 0 S candidate events is shown in Fig. 4. four-pion background (hatched, modeled as a multi-step function). The requirement | cos θ * | < 0.8 is applied.

III. BACKGROUND SUBTRACTION
We first consider non-exclusive background of the type K 0 S K 0 S X, where X is one or more particles. Then we discuss four-track events: π + π − π + π − and K 0 S K ± π ∓ .

A. Non-exclusive background
The contamination by the non-exclusive background process, K 0 S K 0 S X, is estimated by fitting the p t -balance (| p * t |) distribution with a function in which both the signal and background are considered in the region below 0.18 GeV/c. The region above 0.18 GeV/c is not used in this estimate because the p t -balance requirement effectively suppresses events in this region. We approximate the signal distribution with a function that is determined empirically from a signal MC study: where x ≡ | p * t |, α = 1.56 is determined from signal MC, and the parameters A, B and C are floated in the fits in each bin of W and | cos θ * |.
The background distribution is approximated with first-and second-order polynomials connected smoothly at x = 0.05 GeV/c: We verify this approximation in our analyses of the π 0 π 0 and ηπ 0 two-photon production where we observed a large amount of non-exclusive background of the same type [8][9][10]. The fit is performed for data in two-dimensional (W , | cos θ * |) bins of width ∆W = 0.1 GeV (0.2 GeV) for W below (above) 2.0 GeV and ∆| cos θ * | = 0.2.
The results of several such fits are shown for the 0.2 < | cos θ * | < 0.4 region in Fig. 5. The background component is small in the signal region where the data are well described by our parameterization.
To extract the signal yields from data, we subtract the background contributions from our fits. The (W , | cos θ * |) dependence of the background is approximated with a continuous function that is quadratic in most of the W range (connected to linear in a subset of this range) and linear in | cos θ * |. The background yields in each W region, integrated over the angular bins, are shown in Fig. 4.
We estimate the systematic uncertainty associated with this background and its subtraction as half of the subtracted component. We add another 2% error in quadrature to account for the uncertainty in the background p t fit procedure.
Background from the four-pion process is estimated using the summed yield in the M K sideband regions, 0.4826-0.4876 GeV/c 2 and 0.5076-0.5126 GeV/c 2 ; the sum of the widths is the same as that for the signal region (0.4926-0.5026 GeV/c 2 ). We show M K distributions for data in some W regions in Fig. 6. The background contribution is appreciable in the region W < 2.2 GeV only; as this background is always less than 1% for W > 2.2 GeV, we incorporate the uncertainty in our estimate of this contribution in the systematic error but perform no subtraction in this W region.
We obtain the W distribution of the M Ksideband yields for the four separate | cos θ * | bins with a bin width of 0.2. To subtract the four-pion background, we approximate the (W , | cos θ * |) dependence of the background with a multi-step function for W (as shown in Fig. 4) and a linear function for | cos θ * |.
If there were an overlap in the two kinds of backgrounds, i.e., if non-exclusive four-pion events (π + π − π + π − X) were to mimic the K 0 S K 0 S X background, these contributions would be doubly counted and over-subtracted. We find no significantly large non-K 0 S K 0 S (X) contribution in the M K distribution for the p t -unbalanced events with 0.1 GeV/c < | p * t | < 0.2 GeV/c, and therefore estimate the systematic uncertainty associated with the background subtraction as a half of the subtracted component. The possible effect of the overlap is included in this systematic uncertainty.
The K 0 S K ∓ π ± two-photon production, which has a cross section about ten times larger than that of the signal, would contaminate the signal sample if the charged kaon were misidentified as a pion.
According to our MC-based studies, the probability that a generated K 0 S Kπ event is selected as a K 0 S K 0 S signal candidate is smaller than ∼ 10 −4 for W > 2.0 GeV. This probability is so small because of the requirement on the decay vertex distances r V 1 and r V 2 imposed to reject this background.
We use two data-based methods to estimate the remaining K 0 S Kπ background: from a study of the r V distributions near the IP and using our previous measurement of the K 0 S Kπ production process [41]. In the first method, we investigate the r V distribution after identifying one K 0 S with a large r V j , r V j > 1 cm on the opposite side. An excess of events near r V = 0 cm is observed in data for W < 2.0 GeV/c. This is due to the K 0 S Kπ background process, constituting between 0.1% and 4% of the sample at larger r V . This component is observed primarily in the W region below 1.5 GeV. The concentration of the background in the W region may be partially due to four-pion final processes, where one pion track is misreconstructed, resulting in a fake reconstructed vertex. Since we do not separate the four-pion and K 0 S K ∓ π ± backgrounds clearly in the low-W region, we subtract this background assuming the contribution to be 2% ± 2% of the signal in the W region below 1.5 GeV. For W > 1.5 GeV, the excess in the r V distribution is small; this is supported by a study using the measurement of K 0 S Kπ production.
In the second method, the observed yield from the process γγ → K 0 S Kπ is an order of magnitude larger than that of the signal process for W > 2.5 GeV [41], but this background is suppressed by a factor of ∼ 1000 in the data sample after our selection criteria are applied. Thus, it contributes less than 1% to the signal sample. We take this possible contamination into account as a systematic uncertainty of 1% for W > 1.5 GeV.

IV. EFFICIENCY AND EFFICIENCY CORRECTIONS
In this section, we describe efficiency estimates including the factors from the L4 filter, triggers and K 0 S K 0 S reconstruction. Then we discuss corrections for beam energy dependence.
A. The L4 efficiency Some loss of efficiency is introduced by the L4 software filter that is designed to suppress beam-gas and beam-wall events. Figure 7 shows the dependence of the L4 efficiency on W for signal MC events that pass the trigger and all the selection criteria for an assumed isotropic angular dependence. The efficiency is significantly reduced for W < 1.1 GeV and is stable, in the range between 80% and 94%, for W > 1.1 GeV. For very low W , the inefficiency is dominated by the low reconstruction efficiency for tracks with small p t ; for high W , it is explained by tracks with large dr.

Tuning of the simulator for trigger B
We tune the energy threshold for the ECL trigger (LowE), whose nominal value is 0.5 GeV, by comparing the efficiency curves of trigger B between the data and MC events in the four-pion process. With this tuning study in the trigger simulator (TSIM), the optimal value is determined to be 0.52 GeV.
In addition, we find a disagreement of about 20% between data and MC for the energy deposition in the ECL by a low-energy pion. As it is impractical to make dedicated changes in the trigger or detector simulation to describe the detector response to low-energy charged pions for this analysis, we have effectively shifted the LowE threshold by +110 MeV (to 0.63 GeV) to compensate for the pion-energy deposition mismatch.
This shift could affect the efficiency of the selection criterion based on E ECL . We study this possible effect and conclude that it is small, because of the loose criterion on E ECL . As our studies indicate that we could underestimate the efficiency by ∼ 1% because of the ECL energy shift, we correct its value by this amount and assign 1% to the systematic uncertainty of this selection in the entire kinematic region.

Estimation of the trigger efficiency
Using TSIM, we estimate the trigger efficiency for the combination of triggers A, B and C. Its validation using data is non-trivial, because we do not have mutually exclusive triggers to precisely measure the trigger efficiency from the data alone. We find that the contribution of trigger C to the combined efficiency is very small (0.3% -2.0%, depending on W and | cos θ * |), so its contribution to the systematic error is negligible. To estimate the systematic uncertainty of the combined trigger efficiency, we study "trigger- The difference in the angular distribution between the MC and data could cause an apparent deviation of the trigger efficiencies and their ratios in the comparison: in MC, we implement a flat distribution while, in data, steep changes of the distribution are seen for the small angles (typically, in 0.5 < | cos θ * | < 0.8) in some energy regions. To reduce this artifact in the plot for the region 0.4 < | cos θ * | < 0.8 ( Fig. 8(b,d,f)), we subdivide the region into two bins with the same width, 0.2, and take an average for the two bins, for the trigger-A and -B efficiencies and the ratio. The trigger efficiencies estimated by the data and the MC simulation agree within 0.05 except for a low-statistics region.
The assumption of flat angular distributions in MC introduces no bias in the efficiency calculation for cross section derivation because the efficiency is estimated on a bin-by-bin basis with a further narrow bin width, 0.05, in | cos θ * |, whose resolution is much finer than the bin width, as described in Sec. IV E.
In Fig. 7, we show the TSIM trigger efficiency as a function of W for isotropically simulated MC events that satisfy the L4 and all other selection criteria in our analysis. The trigger efficiency rises steeply from 3% near W = 1.05 GeV to 90% near W = 1.6 GeV.
The systematic uncertainty of the trigger efficiency is estimated using the differences in the trigger-A and -B efficiencies and ratios between data and MC, taking into account the correlation between the triggers A and B as estimated from MC. It is evaluated to be 5%-7%, with a weak W dependence.

Validation of the trigger efficiency
We compare our data with the results from the L3 experiment for the cross section of the γγ → 4π process [43], where the π + π − π + π − final-state includes ρ 0 ρ 0 production but not K 0 S K 0 S production. Ideally we would prefer to compare our results directly with K 0 S K 0 S data obtained in previous experiments; however, no such high-statistics data are available. Figure 9 shows a comparison between Belle and L3 for the cross section of the four-pion process (excluding K 0 S K 0 S ) at seven W points (the W bin widths being different between Belle and L3). The Belle selection for the four-pion process has a p t balance cut at 50 MeV/c and non-exclusive backgrounds are subtracted using the p t distribution.
The relative systematic error of the Belle result is estimated to be 10%, while the statistical error is much smaller than that of L3. The Belle result is consistent with the L3 results, but no accurate comparison at a level better than 10% is possible. We assume the L3-determined fractions of the ρ 0 ρ 0 components with spin 0 and 2. Note that the efficiency of the four-pion final state depends on this assumption.
The systematic error associated with the selection efficiency of the K 0 S pairs is estimated by varying the selection criteria in the signal MC. When we do not find two K 0 S candidates with our nominal criteria, we loosen the |∆z| criterion to |∆z| < 10 cm, remove the requirements on K 0 S vertices and loosen the requirement on M K to | M K − m K | < 10 MeV/c 2 , keeping all other criteria at their nominal values. These changes increase both signal efficiency and backgrounds, and we evaluate them with the same methods. The increase of the efficiency is 3%-10% (10%-20%) for W > 1.15 GeV (W < 1.15 GeV).
After the background subtraction, we use the differences in the fractional increase of the efficiency between the original and the loose cuts as its systematic uncertainty. It is difficult to evaluate backgrounds below W < 1.3 GeV because the contamination is larger than the efficiency change and the two different types of non-K 0 S K 0 S backgrounds are not well separated. As the systematic uncertainty is not expected to strongly depend on W , we assign 3% for W < 2.6 GeV and 5% for W > 2.6 GeV as the uncertainty in the efficiency reconstruction for the K 0 S pairs. Figure 10 shows the distribution of cos θ (cosine of the laboratory angle) of K 0 S for the signal candidates at W =1.7 -1.9 GeV for the data and MC. Good agreement between the data and MC is obtained except for the forward-most bin (cos θ > 0.9). The discrepancy there is due to the inadequate efficiency estimation, but its effect (about 3%) is within the systematic uncertainty from tracking, K 0 S reconstruction and trigger efficiencies (see Sec. V).

D. Beam energy dependence
The beam-energy dependence of the luminosity function and the efficiency is studied at the three en- generated at each energy. We compare the luminosity function and efficiency at several W points among the three beam energies. We use 10.58 GeV as the reference energy point and apply a correction proportional to the integrated luminosity to each sample at the other energies.
The luminosity function has a beam-energy dependence with a factor depending on W ; for W in (1.1 GeV, 2.0 GeV, 4.0 GeV), the factor is (−5%, −6%, −10%) for 10.02 GeV and (+2%, +3%, +5%) for 10.88 GeV. Meanwhile, the efficiency depends on the beam energy: +3% at 10.02 GeV and −1% at 10.88 GeV, which is opposite to the trend in the luminosity function. It is also weakly dependent on W .
The overall effect of the beam-energy differences is negligible when we apply the values of the efficiency and luminosity function for 10.58 GeV to all the data, and it is estimated to be at most 0.4% at any W . We do not correct for this effect and do not assign any systematic error.

E. Invariant-mass and angular resolution
We estimate a K 0 S K 0 S mass resolution (i.e., a W resolution) of σ W /W = 0.2% for the entire W region, with a small W dependence, according to a signal MC study. As this is much smaller than the bin width (at worst, σ M < 4 MeV near W = 1.9 GeV, where the bin width is 10 MeV), we do not apply unfolding. The estimated systematic shift due to bin migrations associated with resolution is less than 1% and is absorbed in the systematics.
The resolution for the c.m. angle measurement in each event is typically σ | cos θ * | = 0.0025, which is much smaller than the bin width of 0.1. The measured cross section for the γγ → π + π − π + π − process including ρ 0 ρ 0 from Belle (closed circles) and L3 (diamonds) [43]. The error bars include both statistical and systematic uncertainties, with a uniform 10% estimate used for Belle. These distributions are used solely for trigger-efficiency validation.
candidate events at W = 1.7 -1.9 GeV and | cos θ * | < 0.8 (two entries per event) for the data (dots) and MC (histogram). MC distribution is normalized to have the same number of kaons as observed in data.

V. DIFFERENTIAL CROSS SECTION
The differential cross section dσ/d| cos θ * | is derived after the subtraction of the backgrounds and the application of the corrections to the yields and efficiencies: where N (N bkg ) is the number of candidate (background) events, Ldt is the total integrated luminosity and L γγ is the two-photon luminosity function, calculated as a function of W . ∆W and ∆| cos θ * | are the bin widths, and is the efficiency that includes all trigger/selection effects. The W and | cos θ * | dependence of the overall efficiency is shown in Fig. 11. The efficiency is smaller than 0.14 everywhere in the measurement range. A major cause of the overall efficiency loss is associated with a Lorentz boost of the two-photon system which results in at least one K 0 S falling outside of the detector's acceptance typically more than half of the time. Note that this efficiency loss strongly depends on W and | cos θ * |. We extract the differential cross section in the range | cos θ * | < 0.8 and 1.1 GeV< W < 3.3 GeV, with a W bin width of 10 MeV for W = 1.1−1.9 GeV, 20 MeV for 1.9 -2.4 GeV, 40 MeV for 2.4 -2.6 GeV, and 100 MeV for 2.6 -3.3 GeV. In this extraction, we first evaluate the differential cross section for finer bin widths, ∆W = 10 MeV and ∆| cos θ * | = 0.05 over the entire region, using the efficiency for the central point of each bin. The values for these fine bins are then combined via a weighted average into the coarser bins, with a weight calculated from the statistical errors.
In the range W = 1.05 − 1.10 GeV, we extract only the cross section integrated over | cos θ * | < 0.6, assuming a flat angular dependence of the differential cross section because of limited statistics and the limited coverage in the forward angles in the vicinity of | cos θ * | ∼ 0.6. In the range W = 3.3−3.6 GeV, we do not extract the γγ → K 0 S K 0 S cross section where the contributions from the χ c0 and χ c2 resonances dominate the yield; we cannot subtract leakages from these narrow states reliably over the entire region.
In the range W = 3.6 − 4.0 GeV, we find some contribution from the signal process. It is possible to extract the integrated cross section for | cos θ * | < 0.8 in this W region; however, we do not present differential cross section due to small statistics. There could be a contribution from high-mass charmonium resonance(s) (χ cJ (2P ) for example) at W = 3.80 − 3.95 GeV, as we find some events at large angles in this W range; these events are included in the total cross section (see Fig. 28).
At W = 4.0−4.1 GeV, we find only a small number of signal events that give a peak near | p * t | = 0, consistent with a large background contamination. No cross section measurement is therefore performed in the W region above 4.0 GeV. Figure 12 shows the cross section integrated over | cos θ * |. The integration is performed by summing the differential cross section for | cos θ * | < 0.8 or | cos θ * | < 0.6. The error bars are statistical only. The curves in the figure show the total systematic errors.
The systematic error includes contributions from the uncertainties in tracking efficiency (2% for 4 tracks), beam-background effects (1%) estimated from the stability of yield ratios between the data and MC across all run periods, pion identification (2% for four pions), non-exclusive and four-pion backgrounds (described in Sec. III A and B), geometrical coverage and fit uncertainty (4% in total), K 0 S Kπ background subtraction (Sec. III C), K 0 S -pair reconstruction (Sec. IV C), trigger efficiency (Sec. IV B), and the E ECL cut (Sec. IV B). We assign the uncertainty for the L4 efficiency to be about 10% of the inefficiency in different W regions. The systematic error associated with the uncertainty in the integrated luminosity and luminosity function includes the form-factor effect of space-like photons. Summing in quadrature, the total systematic uncertainty evaluated is typically 10%. The systematic uncertainties are summarized in Table I. Figure 12(a) shows the measured integrated cross section (| cos θ * | ≤ 0.8), where prominent peaks are observed near 1.3, 1.5 and 1.8 GeV. Enhancements are also observed around 2.3 and 2.6 GeV. A close-up view of the integrated cross section (| cos θ * | ≤ 0.6) near the threshold is shown in Fig. 13, where the cross section is small (< 1 nb), in agreement with theoretical predictions [15,16].

VI. STUDY OF RESONANCES
In this section we describe the extraction of partial wave information from our data by fitting the differential cross section using suitable parameterizations to estimate the mass, total width and Γ γγ B(KK) of the f 2 (1525), to derive the phase difference between the f 2 (1270) and a 2 (1320) and to identify the nature and obtain the parameters of the resonances near 1.8, 2.3 and 2.6 GeV.

A. Partial wave amplitudes
In the γγ → K 0 S K 0 S channel, only the partial waves of even angular momenta contribute. Furthermore, in the energy region W < ∼ 3 GeV, the J > 4 partial waves may be ignored, so only S, D and G waves are considered in the fit. The differential cross section can then be expressed as where S is the S-wave amplitude, D 0 and G 0 (D 2 and G 2 ) denote the helicity-0 (2) components of the D and G wave [44], respectively, and Y m J are the spherical harmonics. The angular dependence of the cross section is governed by the spherical harmonics while the energy dependence is determined by the partial waves.
Since the spherical harmonics are not independent of each other, a unique decomposition of partial waves cannot be determined using the measured differential cross section. To overcome this problem, we rewrite Eq. (7) as The "hat-amplitudes"Ŝ 2 ,D 2 0 ,D 2 2 ,Ĝ 2 0 andĜ 2 2 can be negative because of interference terms, and correspond to |S| 2 , |D 0 | 2 , |D 2 | 2 , |G 0 | 2 and |G 2 | 2 , respectively, when interference terms are ignored [8]. The W dependence of the γγ → K 0 S K 0 S cross section after integrating over the angle up to (a,b) | cos θ * | < 0.8 (black points) and (c,d) | cos θ * | < 0.6 (blue points). The orange square markers in (d) are from our previous publication [5] for | cos θ * | < 0.6. The solid curves are the systematic uncertainties.
As the absolute squares of the spherical harmonics are independent of each other, we can fit the differential cross section in each W bin to obtainŜ 2 ,D 2 0 , D 2 2 ,Ĝ 2 0 andĜ 2 2 . The fit with the value of J ≤ 4 is named the "SDG fit." At low energy, we expect that the contribution from J = 4 is negligible, so we also perform a separate fit by settingĜ 2 0 =Ĝ 2 2 = 0, which is named the "SD fit." The differential cross section is fitted according to Eq. (8), where statistical errors only are taken into account. The differential cross section for | cos θ * | ≤ 0.8 is extracted for 1.1 GeV ≤ W ≤ 3.3 GeV. In the SDG fit, two consecutive data points of ∆W = 0.01 GeV are merged, resulting in bins of width 0.02 GeV.
The obtained spectra ofŜ 2 ,D 2 0 andD 2 2 for the SD fit are shown in Fig. 14. Figures 15 and 16 show the hat-amplitudes for the SDG fit.Ĝ 2 0 ±Ĝ 2 2 are also plotted in Fig. 16, since the angular dependences of |Y 0 4 | 2 and |Y 2 4 | 2 are similar for | cos θ * | < ∼ 0.6. In the SDG fit, the structures inD 2 2 are less visible and the G waves appear to be small for W ≤ 3.3 GeV. Figure 17 shows the differential cross section for selected W bins with the fittedŜ,D 0 andD 2 waves.
Although the derived hat-amplitudesŜ 2 ,D 2 0 ,D 2 2 , G 2 0 andĜ 2 2 in fact contain interference terms such as (S * D 0 ), they do provide useful information about partial wave amplitudes. Two prominent peaks are observed in theD 2 2 spectrum: the peak near 1.3 GeV  is due to the interference between the f 2 (1270) and a 2 (1320) and the second peak is due to the f 2 (1525). No other notable structures are observed in Fig. 14 (right). In theŜ 2 spectrum shown in Fig. 14 (left  top), three peaks around 1.8, 2.3 and 2.6 GeV are observed. The lowest may be due to the f 0 (1710) (not a tensor meson, as discussed in past experiments [4,17,21]). This might be an a 0 meson, though no such mesons have been observed previously in this mass region [23].D 2 0 is rather small and featureless except around 2.1 and 2.6 GeV, and hence the D 0 wave may also be small but non-zero: there appears to be an interference term between S and D 0 .
We use our assumptions for the partial wave amplitudes and fit the data to extract the parameters of the resonances. Note that we do not fit the obtained spectra of hat-amplitudesŜ 2 ,D 2 0 andD 2 2 , but rather fit the differential cross section directly using Eq. (7). In our analysis, we fit the energy region W ≤ 3.0 GeV, with separate fits for W ≤ 2 GeV and W > 2 GeV. In this section, we describe our fits in the W ≤ 2.0 GeV region. Motivated by the spectra ofD 2 2 and S 2 , we include the f 2 (1270), a 2 (1320) and f 2 (1525) in the D 2 wave and test the hypothesis of a scalar meson (coined the f 0 (1710)) in the S wave. In this analysis, we measure the relative phase probing the destructive interference between the f 2 (1270) and a 2 (1320) and determine relevant parameters of the f 2 (1525), in particular, Γ γγ B(KK).

Parameterization of amplitudes
Based on the above observation, the amplitudes S, D 0 and D 2 are parameterized as follows: where A f0 (1710 φ D0 and φ D2 are the phases of the resonances and background amplitudes. The phases are defined relative to B S (f 2 (1270)) for helicity-0 (2) amplitudes. Using this convention, the relative phase between the f 2 (1270) and a 2 (1320) is given by φ a2 . We also study the case in which the f 0 (1710) is replaced by a tensor meson (labeled the f 2 (1710) here, although the f 2 (1810) is listed in PDG [23]) in D 2 . To investigate if our approximation could describe the data well without this resonant contribution, we also perform a fit assuming no resonance at 1.8 GeV.
We assume the background amplitudes to be quadratic in W multiplied by the threshold factor β 2l+1 for all waves: the velocity of the K 0 S divided by the speed of light, m K 0 S is the mass of the K 0 S , and l is 0 (2) for S (D 0 and D 2 ).
We use the parameterization of the f 2 (1270) and f 2 (1525) given in Ref. [8] and that of the a 2 (1320) in Ref. [10]. We note that B(R → K 0 The amplitude A R (W ) for each spin-J resonance R of mass m R is parameterized using the relativistic

Breit-Wigner formula
Hereinafter, we consider the case J = 2. The energydependent total width Γ tot (W ) is given by where X 1 , X 2 is π, K, η, γ, etc. For J = 2, the partial width Γ X1X2 (W ) is parameterized as [45]: where Γ R is the total width at the resonance mass, and r R is an effective interaction radius that varies from 1 (GeV/c) −1 to 7 (GeV/c) −1 in different hadronic reactions [46][47][48]. For the three-body and other decay modes,  is used instead of Eq. (13). This formalism is used for the f 2 (1270), a 2 (1320) and f 2 (1525). All parameters of the f 2 (1270) and a 2 (1320) are fixed at the PDG values [23] except for r R , which is fixed at the value determined in Refs. [1,2], as summarized in Tables II  and III. Finally, the parameterization of the f 0 (M ) meson for M = 1710 MeV/c 2 is taken to be: where Γ γγ B(KK) f0 is the product of the two-photon width and the branching fraction to KK for the f 0 (M ) meson. Its PDG [23] parameters are summarized in Table IV, together with the parameters (when known) of the f 2 (1810) and a 2 (1700).

Fit in the region 1.2 GeV ≤ W ≤ 2.0 GeV
We perform a fit for the region 1.2 GeV ≤ W ≤ 2.0 GeV by floating the mass, width, Γ γγ B(KK) and  Twenty parameters describing the assumed amplitudes are obtained by fitting the differential cross sections. To search for the global minimum goodness of fit χ 2 min to identify possible multiple solutions, about 1000 sets of randomly generated initial parameters are employed for fits performed using MI-NUIT [49]. A fit is accepted as a satisfactory solution if its χ 2 -value is within χ 2 min + 10 (corresponding to 3.2σ).
If the f 0 (1710) hypothesis is assumed to explain the peak at W ∼ 1.8 GeV, four solutions are obtained with χ 2 min /ndf = 677.2/580, where ndf is the number of degrees of freedom. These solutions are distinguished by the Γ γγ B(KK) value, which ranges from 6.3 to 216 eV for the f 0 (1710), and from 83 to 104 eV for the f 2 (1525), as listed in Table V.
When the f 2 (1710) hypothesis is used, the two obtained solutions have lower quality with χ 2 min /ndf = 755.6/580. Their fitted values are also listed in Table V. As the f 0 (1710) solutions have lower χ 2 min /ndf , they are favored over the f 2 (1710). We conclude that the region 1.2 GeV ≤ W ≤ 2.0 GeV is too wide to fit in extracting the desired parameters at once. We therefore fit individual parameters one at a time, keeping in mind the limitations of this approach. 3. The "f 2 (1525) fit" Based on the above observation, we first obtain the f 2 (1525) parameters by fitting the c.m. energy range 1.15 GeV ≤ W ≤ 1.65 GeV and ignoring the contribution of the f J (1710). The differential cross section is fit with the parameterized amplitudes by floating the f 2 (1525) parameters as well as the rela- tive phase between the a 2 (1320) and f 2 (1270). Hereinafter, this fit is referred to as the "f 2 (1525) fit." The background amplitudes are approximated with linear functions because the fitting range is rather narrow. There are thirteen parameters to extract from this fit. Two solutions are obtained, both with χ 2 /ndf = 0.97. The main difference between the two solutions is the values of Γ γγ B(KK) for the f 2 (1525): 113 and 48 eV, with the two solutions referred to as H (high) and L (low), respectively. They correspond to destructive and constructive interference between the f 2 (1525) and non-resonant D 2 background. The fitted results are shown in Figs. 18 and 19 for the H and L solutions, respectively. The fitted values are listed and compared with those of PDG [23] in Table VI. The quoted errors are MINOS statistical errors, determined by evaluating the parameter values that give χ 2 min +1 for each variable being studied. In the fit, all other parameters are floated. In both solutions, the interference between the f 2 (1270) and a 2 (1320) is indeed destructive as predicted [13], i.e., the measured φ a2 is close to 180 • .
The following sources of systematic uncertainties for the fitted parameters are considered: dependence on the fit region, normalization errors of the differential cross section, assumptions on the background amplitudes and assumed parameters of the f 2 (1270) and a 2 (1320). In each study, a fit is performed that allows all the parameters to float; the differences of the fitted parameters from the nominal values are quoted as systematic uncertainties for both solutions, H and L.
Two fitting regions shifted by ±0.05 GeV (10% of the W -range), are used to estimate the systematics associated with the fitting range. The systematic uncertainties associated with the normalization are separated into two groups: one from the overall normalization (±4.0%) and the other from the distortion of the spectra in both | cos θ * | and W . To estimate the uncertainty associated with the overall normalization, fits are performed by shifting the cross sections coherently by ±4%. For point-by-point normalization, fits are performed by shifting the cross section by ±|dσ/dΩ| × σ (W,| cos θ * |) , where σ is the relative uncertainty of the efficiency (referred to as Efficiency in Table VII). For the spectral distortion studies, the differential cross sections are shifted by ±0.1 × |dσ/dΩ| × (| cos θ * | − 0.4)) (referred to as | cos θ * | bias) and ±0.08 × |dσ/dΩ| × (W − 1.4 GeV) (referred to as W bias). We use the absolute value of dσ/dΩ because some of the central values for measured differential cross sections are negative due to background subtraction.
For studies of background (BG) amplitudes, each background wave is approximated by a constant or a parabola. The value of r R is changed by ±0.03 (GeV/c) −1 according to Refs. [1,2]. Finally, the parameters of the f 2 (1270) and a 2 (1320) are changed one by one by their uncertainties shown in PDG [23].  The total systematic uncertainties are calculated by adding the individual uncertainties in quadrature. The resulting systematic uncertainties are summarized in Table VII. In some of our studies, the value of Γ γγ B(KK) for the f 2 (1525) fluctuates between the H and L solutions. The obtained results for the relative phase between the a 2 (1320) and f 2 (1270) and parameters of the f 2 (1525) are given in Table VI.

A fit including the fJ (1710)
We fit the region 1.2 GeV ≤ W ≤ 2.0 GeV by fixing the parameters of the f 2 (1525) and φ a2 to either the H or L solution, and by including the contribution of the f 0 (1710) (coined the "f 0 (1710) fit"). The background amplitude is assumed to be a secondorder polynomial, whose parameters are floated in the fit.
A unique solution is obtained for each of the H and L solutions (named "fit-H" and "fit-L," where H and L stand for the H and L solutions of the f 2 (1525) fit, respectively). These solutions are summarized in Table VIII. Figures 20 and 21 show the fitted results for fit-H and fit-L, respectively. Figures 22 and 23 show fit-H and fit-L solutions superimposed on the differential cross section for selected W bins.
We also study a case where the structure near W = 1.8 GeV is assumed to be due to a tensor meson (labeled the f 2 (1710), which can be either a 2 (1700) or f 2 (1810) (referred to as the "f 2 (1710) fit"). The contribution from tensor mesons may be suppressed due to destructive interference between the f 2 (1810) and a 2 (1700); this hypothesis could also be tested by analyzing γγ → K + K − data. A unique best fit with poor χ 2 is obtained for the f 2 (1710) fit with either of the H and L solutions of the f 2 (1525) fit. Thus, the hypothesis of J = 2 for the f J (1710) is disfavored by the data. Fitted values are summarized in Table VIII. Figures 24 and 25 show the fitted results for the f 2 (1710) fit for each of the H and L solutions.
Furthermore, we fit the hypothesis where we assume no resonance near W = 1. hypothesis with even worse χ 2 /ndf of 1349.8/589. We conclude that our fit favors the presence of the f 0 (1710) in our data.
Systematic uncertainties are estimated similarly to those for the f 2 (1525) fit. In the W -range study, fits are performed in two fit regions: 1.12 GeV ≤ W ≤ 1.92 GeV and 1.28 GeV ≤ W ≤ 2.08 GeV. For W -distortion, a study is performed by shifting the cross section by ±0.08 × |dσ/dΩ|(W − 1.6 GeV); for background waves, by changing each wave to a firstor third-order polynomial; and for the parameters of the f 2 (1525), by shifting the values by their MINOS errors. The results for the systematic uncertainties are summarized in Table IX. Table VIII lists the results for the f 0 (1710) fit (fit-H and fit-L).

Final results for the region W < 2.0 GeV
As described above, we obtain two solutions referred to as H and L for the f 2 (1525) fit, and corresponding fits are performed in the f 0 (1710) fit by fixing the f 2 (1525) parameters to those of either the H or L solution. Here, we combine solutions statistically to obtain final results.
From each pair of solutions for an observable x, a probability density function (PDF) P (x) is formed as the sum of asymmetric Gaussian functions that correspond to the two solutions with asymmetric errors. These functions are weighted according to the χ 2 differences between the two solutions. The most probable value x f is the one that gives the maximum in P (x). Asymmetric statistical errors σ l and The systematic uncertainty for observable x is determined from the solution with the largest deviation from x f . The final results thus obtained are listed in Tables VI and VIII. C. Fitting the region W > 2.0 GeV We investigate the structures around 2.3 and 2.6 GeV. In fitting the region 2.0 GeV ≤ W ≤ 3.0 GeV, we assume that the non-resonant backgrounds in the S, D 0 , D 2 , G 0 and G 2 waves obey a power law in W . When we parameterize them using a polynomial approximation as in Eq. (10), we obtain fits of poor quality. We parameterize the backgrounds as where the index i denotes S, D 0 , D 2 , G 0 or G 2 waves, W 0 is chosen to be the lower boundary of the fitting region (nominally W 0 = 2.0 GeV), and b i and c i are the free parameters. The phases of B S and B D2 are chosen to be zero as a reference for the other phases. The parameters b i are set positive to resolve arbitrary sign ambiguities. We also investigate a possible contribution from the J = 4 resonances. Table X summarizes the parameters of the f 0 (2200), f 2 (2300) and f 4 (2300) that are known to couple to KK [23]. We allow B G0 and/or B G2 to be non-zero in Eq. (18). We fit 13 assumptions for the structures around 2.3 and 2.6 GeV that are observed in the plot of the integrated cross section shown in Fig. 12(a). A fit performed assuming the presence of the f J (2200) (J = 0, 2, 4) and f J (2500) (J = 0, 2, 4) is referred to as an "f J -f J fit." We also investigate hypotheses in which there are no resonances (or only one) for the two structures. These fits are referred to as "no-resonance-" ("only-f J -") fits, respectively.
When both B G0 and B G2 are allowed to be nonzero, too many solutions are obtained because of the several combinations of interfering amplitudes (not shown). Thus, we focus on the hypothesis wherein only B G2 is non-zero. This choice is based on the idea that the possible resonances f 4 (2200) and f 4 (2500) are included in the G 2 wave only because of helicity considerations. A summary of the fitted results is given in Table XI. In this case, once again, some of the f J -f J fits give multiple solutions. In some cases, one or more c i values in Eq. (18) assume unphysically large values. Thus, we constrain the maximum values of c i to be 20.
A unique solution of relatively good quality is obtained for the f 2 -f 0 fit, while other hypotheses yield larger values of χ 2 /ndf . The f 2 -f 0 fit is also favored for the case in which both B G0 and B G2 are assumed to be non-zero. Thus, we conclude that the structure around 2.3 GeV is likely due to a tensor meson (referred to tentatively as f 2 (2200)) and the one near 2.6 GeV is likely to be a scalar meson (possibly f 0 (2500)).
The fitted values obtained from the f 2 -f 0 fit are summarized in Table XII for the mass, total width     In this subsection, we summarize the fitted results and discuss their implications. First, the destructive interference between the f 2 (1270) and a 2 (1320) [13] is confirmed with high accuracy; the relative phase, φ a2 , is measured to be 172.6 +6.0+12.

+108
−12 eV, respectively. The systematic uncertainty of Γ γγ B(KK) is fairly large. Nevertheless, this is the first measurement of this parameter that includes the interference with a non-resonant amplitude.
The structure near 1.6 GeV is attributed to a scalar meson and is interpreted to be the f 0 (1710). To obtain the significance for the assignment of the f 0 (1710) over that of the f 2 (1710), fits are performed for each of the sources of systematic uncertainty for the two hypotheses and the minimum χ 2 difference is identified among these fits. It is found to be 63.3, which corresponds to a significance of 7.9σ favoring the f 0 (1710).
A similar study is performed for the f J (2200) hypothesis by comparing the values of χ 2 min of the f 2 -f 0 fit and f 0 -f 0 fit for each source of systematic uncertainty. We obtain a minimum χ 2 difference to be of 11.3, corresponding to a 3.4σ significance. For the f 0 (2500), the f 2 -f 0 fit gives the best χ 2 ; the next best, the f 2 -f 2 fit, yields a 4.3σ significance. Thus, while we cannot make definitive assignments about the spins of the f J (2200) and f J (2500), J = 2 and J = 0 hypotheses, respectively, are favored.
The values of Γ γγ B(KK) for the f 0 (1710), f 2 (2200) and f 0 (2500) are estimated for the first time and are found to be 12 +3+227 −7−40 eV, respectively. Each value could provide important information on the constituents of the corresponding resonance. For example, the f 0 (1710) is identified as an unmixed scalar glueball according to a coupled-multi-channel analysis [50]. However, the f 0 (1710) is unlikely to be a glueball candidate because the observed value of Γ γγ B(KK), combined with the implied value of Γ γγ B(ππ) ( Γ γγ B(KK) for the flavor-SU(3)-symmetric decay of a glueball) would indicate a large two-photon width, contrary to the expectation of much less than 1 eV for a glueball (see, e.g., Refs. [26][27][28][29]). Therefore, we conclude that the f 0 (1710) is a resonance with a large ss admixture.
The measured mass of the f 2 (2200), 2243 +7+3 −6−29 MeV/c 2 , is close to that of the f J (2220) and smaller than that of the f 2 (2300) [23]. The f J (2220) is usually assumed to be a glueball candidate due to the small value of its total width (23 +8 −7 MeV). The structure found by Belle in the γγ → K + K − reaction around 2.3 GeV [4] is interpreted as the f 2 (2300) by PDG [23]. The measured total width of the f 2 (2200), 145 ± 12 +27 −34 MeV, is much wider than that of the f J (2220) and is similar to that of the f 2 (2300).

VII. DERIVATION OF CHARMONIUM CONTRIBUTION
We present our estimates of the χ c0 and χ c2 contributions by measuring the yields of the fit components in the region | cos θ * | < 0.5 and 2.8 GeV < W < 4.0 GeV (Fig. 28). We use only samples with | cos θ * | < 0.5 to enhance the fraction of the charmonium components while disentangling them from the   continuum contribution. We derive Γ γγ (R)B(R → K 0 S K 0 S ) for these charmonium states. We also search for a possible contribution from a higher-mass charmonium in the 3.6 -4.0 GeV/c 2 region.
We assume the angular distributions for the χ c0 and χ c2 components to be flat and proportional to sin 4 θ * (from pure helicity-2 [51]), respectively, to derive Γ γγ B from the yield in | cos θ * | < 0.5. We discuss the effect of interference with the continuum.  A. Evaluation of parameters for the χcJ charmonia The peak structures observed in the yield distribution for 3.3 GeV < W < 3.6 GeV (Fig. 28) are from charmonium production: γγ → χ c0 , χ c2 → K 0 S K 0 S . We fit the distribution to contributions from the χ c0 and χ c2 as well as a smooth continuum component represented by in the W and | cos θ * | ranges 2.9 -3.7 GeV and | cos θ * | < 0.5, respectively, where BW χ cJ (W ) is a Breit-Wigner function for the charmonium amplitude, which is proportional to 1/(W 2 − M 2 χ cJ − iM χ cJ Γ χ cJ ) and is normalized by |BW χ cJ (W )| 2 dW = 1. The component αW −β corresponds to the contribution from the continuum, with a fraction k that interferes with the χ c0 amplitude with a relative phase φ.
It is impossible to determine the interference parameters for the χ c2 by any fits because of its smaller intrinsic width compared to the mass resolution. We fit the χ c2 yield (N χc2 ) with a function in which no interference term is included, as shown by Eq. (19); later, we estimate the maximum effects from the in-  terference term when we evaluate the uncertainty for the two-photon decay width of χ c2 . All parameters except the width of the χ c2 are free. The χ c2 width is fixed to 2.0 MeV, which is smaller than the estimated mass resolution of ∼ 7 MeV. Smearing effects due to a non-zero mass resolution are taken into account in the fit, using a Gaussian function with σ = 7.0 MeV. The small W dependence of the product of the efficiency and luminosity function, ∓0.95% for a change in W of ±10 MeV, is folded in the χ c0 resonance curve.
A binned maximum likelihood method is applied with the bin width ∆W = 10 MeV. We examine two cases: with and without the interference of the χ c0 . We could not determine the k parameter; that is, any 0 < k ≤ 1 gives exactly the same fit quality for the constructive (φ ≈ π/2) and destructive (φ ≈ 3π/2) interference cases. Therefore, the statistical error range for the yield of χ c0 corresponds to the full range of the interference assumption 0 < k ≤ 1.
The maximum effect of the interference of χ c2 with the continuum component is calculated from Eq. (20) because we cannot measure it from the line shape of the charmonium, so we include its maximum possible range in the statistical error. The number of χ c2 events that is proportional to the square of the resonance amplitude is converted from the observed number N χc2 to that with the maximum interference effect N χc2 using the relation where Γ and n I are the total width of the χ c2 and the observed yield density of the continuum component per unit energy in the W range around the χ c2 , respectively.
The fitted results are summarized in Table XIV, where L is the likelihood value. The normalization N χc0 in Eq. (19) is proportional to the square of the resonance amplitude, even when the interference is assumed. The yields from the fits are translated to products of the two-photon decay width and the branching fraction, Γ γγ (χ cJ )B(χ cJ → K 0 S K 0 S ), shown in Table XV.
To estimate the systematic uncertainties associated with the choice of the signal shape approximation, we vary their shape parameters. We change the W resolution from 7 to 8 MeV and modify the term in the denominator of the Breit-Wigner function, from −iM Γ to −iW Γ. The observed changes of the central values of the χ c0 and χ c2 yields are less than 3%. This is because the χ c0 and χ c2 contributions are well separated from each other, and the continuum contribution is very small around the charmonium peaks. The systematic uncertainty is thus dominated by the contributions from the efficiency and luminosity function, and is about 10% in total. The systematic uncertainties for Γ γγ B are shown in Table XV.
The present results are consistent with and supersede those from the previous measurements [5,23]. The interference effect was not taken into account in the previous Belle result.

B. Possible higher-mass charmonium states
We could expect a contribution from the possible higher-mass charmonium states, χ cJ (2P ) (J = 0, 2)  in the W region above 3.8 GeV. The χ c2 (2P ) has been found near 3927 MeV/c 2 in the two-photon process [23] in its decay to the DD state, but the χ c0 (2P ) has not yet been observed in this decay mode. Although no theoretical predictions are available for the branching fractions B(χ cJ (2P ) → K 0 S K 0 S ), a yield of a few events is expected if Γ(χ cJ (2P ) → K 0 S K 0 S ) ≈ Γ(χ cJ (1P ) → K 0 S K 0 S ) and postulated or observed values for Γ tot and Γ γγ for such states are taken.
As seen in Fig. 28, we find 8 events in the W region between 3.7 and 4.0 GeV, consistent with 5.2 events expected in the region from the extrapolated continuum background determined by the fit below 3.7 GeV (with a continuum yield density of dY /dW = 59.2(W/3.5GeV) −13.5 [GeV −1 ], including interference). In the W region between 3.80 and 3.95 GeV, where we expect the presence of contributions from the two higher-mass charmonium states, 7 events are observed, while only 2.3 events are expected from the continuum. The probability for this observation (p-value) is 0.9%.
We evaluate an upper limit for Γ γγ (χ c2 (2P ))B(χ c2 (2P ) → K 0 S K 0 S ). We find 2 events in the χ c2 (2P ) mass region, 3.879 -3.975 GeV/c 2 , which is defined by M ± 2Γ using the known mass and total width [23]. We adopt N U L χc2(2P ) = 5.32 as the upper limit of the yield with a 90% confidence level (CL) for the contribution of the χ c2 (2P ), assuming no background contribution for a conservative limit and based on the Poisson distribution with this mean value giving a 10% probability for two or fewer observed events. This translates into an upper limit for the product of the two-photon decay width and the branching fraction of the χ c2 (2P ) of Γ γγ (χ c2 (2P ))B(χ c2 (2P ) → K 0 S K 0 S ) < 0.064 eV (90% CL) without interference. This upper limit takes into account the uncertainty of the efficiency by increasing the limit by one standard deviation.
C. Search for the decay ηc → K 0 The decay η c → K 0 S K 0 S violates both P and CP invariance. We search for this decay mode in the present data. Copious production of the η c in twophoton collisions has been established in several decay modes by previous measurements [23].
A small peak-like structure near 2.99 GeV seen in Fig. 28 is not statistically significant and corresponds to a fluctuation at the 1.7σ level, which is evaluated from the difference between log-likelihoods for the fits without and with a contribution of the η c , taking into account the interference effect that is described below. We thus set the upper limit of the branching fraction for this decay mode.
We fit the event distribution in the range 2.8 GeV< W < 3.3 GeV with a function similar to Eq. (19) in which the χ c0 contribution is replaced by that of the η c and the χ c2 term is not included. We fix the mass and width of the η c to be 2981 MeV/c 2 and 30 MeV, respectively. The best fit without interference gives N ηc = 5.4 ± 5.0. This is consistent with zero. We determine the 90% CL upper limit with the N UL ηc value that corresponds to the (1.64) 2 worse log-likelihood −2 ln L than that of the best fit.
We take into account uncertainties in the mass, width and the mass resolution associated with our measurement, and repeat the fit by adjusting these parameters by ±2 MeV/c 2 , ±4 MeV and in 5 − −7 MeV, respectively, and choose the most conservative upper limit. The obtained upper limit is N UL ηc = 15 (N UL ηc = 85) without (with) interference. The curves describing the results of the fits used to estimate the upper limits as described are shown in Fig. 29.
The 90% CL upper limits for Γ γγ (η c )B(η c → K 0 S K 0 S ) and B(η c → K 0 S K 0 S ) are summarized in Table XVI; for the latter, Γ γγ (η c ) = 5.3 ± 0.5 keV is used [23]. These upper limits take into account the uncertainties from systematic error of the measurement and the Γ γγ (η c ) value by shifting the limits by a ratio corresponding to 1σ in the direction of increased values.

VIII. QCD STUDIES IN THE HIGH-ENERGY REGION
In this section, the cross-section behavior is studied and compared with predictions from QCD-based models and calculations in the region W > 2.6 GeV. First, we compare the differential cross section with the 1/ sin 4 θ * dependence. Then the W −n behavior of the integrated cross section is examined.
A. Angular dependence of the differential cross section We compare the angular dependence of the differential cross section with the 1/ sin 4 θ * dependence, which is claimed by the handbag model [25]. Earlier Belle measurements for this process supported such a dependence in the W region between 2.4 and 3.3 GeV for | cos θ * | < 0.6 [5].
To make a quantitative statement about the behavior of the cross section, we fit the differential cross section using the approximation A/ sin α θ * , i.e., in each W bin. We summarize the fitted results for the 12 regions in Fig. 30, where the right scales are differential cross sections normalized to the integrated cross section in the range | cos θ * | < 0.8 (that gives the average 1/0.8 = 1.25). This scale is added to improve the visibility of the plots for different W bins. The χ 2 /ndf values obtained from the fits are between 3/6 and 19/6. The obtained W dependence of the parameter α is shown in Fig. 31. The parameter α is found to be above 4 for the W range between 2.7 and 3.3 GeV, but no tendency toward 4 is observed in the high-energy part of the W region. We note that we find a resonance-like contribution considered to be a scalar at around 2.5 GeV, as described in Sec. VI C, which could affect the W dependence of α in the region around 2.4 -2.7 GeV. Information on the meson (M ) distribution amplitude (DA) φ M can be obtained by comparing the observed angular dependence to that of the theoretical calculation [24]; the angular dependence of the data is steeper and more forward-peaked, which indicates that the DA is flatter than assumed.
The function proportional to 1/ sin 4 θ * + b cos 2 θ * that has been applied in our analysis of the γγ → π 0 π 0 process yields fits of poor quality in this study, as the rise of the cos 2 θ * term for the forward angles is insufficient to describe the trend observed in data.

B. W dependence
The W dependence of the cross section integrated over the angle provides important information about the mechanism of the exclusive meson-pair production. We fit the cross section with σ(| cos θ * | < 0.8) = aW −n , for the W region 2.6 -4.0 GeV, excluding 3.3 -3.6 GeV. We exclude the region below 2.6 GeV because a resonance-like contribution is found there. We obtain n = 11.0 ± 0.4. This result is shown by the dashed line in Fig. 32(a). The error is statistical only. We also try the fits for the narrower W region, 2.6 -3.3 GeV, for σ(| cos θ * | < 0.8) and σ(| cos θ * | < 0.6), and obtain n = 10.0 ± 0.5 and n = 11.8 ± 0.6, respectively.
In our previous work, we obtained n = 10.5 ± 0.6 for W = 2.4 − 4.0 GeV excluding 3.3 -3.6 GeV and  The left (right) vertical scale of each subfigure corresponds to the absolute scale (normalized in such a way that the average is 1.25, as described in the text) of the differential cross section.
| cos θ * | < 0.6 [5]. The present analysis in this region yields n = 10.8 ± 0.2. We quote this number only for verification, as we now know that it includes a resonance-like contribution around 2.5 GeV. These results are summarized in Table XVII. We estimate the systematic uncertainty for the n measurement as follows: since the overall normalization error does not affect the determination of n, we consider the W -dependent distortion effect only. As in the resonance studies, we assume ±4% distortions at the two ends of the fit range and continuous variations between them. The distortion for a W range changes the n value with ∆n ≈ log 1.08/ log(W 2 /W 1 ), where W 1 and W 2 de- limit the fit region (chosen to be 2.65 GeV and 3.25 GeV, respectively). We obtain the estimated systematic uncertainty ∆n = 0.4. The slope parameter n that ranges between 10 and 11 for the present process is larger than 6 and 7-8 that are predicted [24] and observed [3], respectively, for the π + π − and K + K − processes. For the process γγ → K 0 S K 0 S , as discussed in Refs. [54,55], the coefficient of the leading-term amplitude is much smaller than that of the non-leading term. Therefore, at this energy the W dependence of the cross section is mainly determined by that of the non-leading terms. In the W region measured in this experiment and including a non-leading term, the perturbative QCD predicts n = 10 [55], which is in reasonable agreement with our measurement.

IX. SUMMARY AND CONCLUSION
We have measured the cross section for the process γγ → K 0 S K 0 S for 1.05 GeV ≤ W ≤ 4.00 GeV with the Belle detector at the asymmetric-energy KEKB collider. The data sample of 972 fb −1 is three orders of magnitude larger than in the previous measurements [17][18][19][20][21]. The differential cross section is measured up to | cos θ * | = 0.8, which allows highsensitivity studies of the amplitudes.
In our study, the differential cross section has been fitted to obtain information on partial waves. The obtained spectra ofŜ 2 ,D 2 0 andD 2 2 indicate the presence of the f 0 (1710), f J (2200) and f J (2500) in addition to the well known f 2 (1270), a 2 (1320) and f 2 (1525). Then fits to the differential cross section are performed by assuming possible resonances in the partial waves. First, we perform a fit in the region 1.15 GeV ≤ W ≤ 1.65 GeV to determine the parameters of the f 2 (1525) as well as the relative phase between the f 2 (1270) and a 2 (1320). Two solutions are obtained and combined statistically. The phase difference between the a 2 (1320) and f 2 (1270) is measured to be 172.6 +6.0+12.2 −0.7−7.0 • , confirming the destructive interference between the two mesons and agreeing with theoretical predictions [13]. The mass, total width and Γ γγ B(KK) of the f 2 (1525) are measured to be 1525.3 +1.2+3.7 −1.4−2.1 MeV/c 2 , 82.9 +2.1+3.3 −2.2−2.0 MeV and 48 +67+108 −8−12 eV, respectively. Note that no interference effect was taken into account in the previous measurements [17][18][19]21].
Evidence for the existence of the f 0 (1710), f 2 (2200) and f 0 (2500) in this channel is obtained. We conclude that the f 0 (1710) and f 2 (2200) are unlikely to be glueballs because their total widths and Γ γγ B(KK) values are much larger than those expected for a pure glueball state.
Analyses in the region W > 2.6 GeV are updated; parameters of the χ c0 and χ c2 and the exponents α and n in (sin θ * ) −α and W −n describing the angular and W behavior of the cross section are extracted from data. The value of α does TABLE XVII: Results for the slope parameter n from the power fit σ ∼ W −n for γγ → K 0 S K 0 S in different fit ranges. The result from the previous work [5] is also shown. The first and second errors are statistical and systematic, respectively. not show the tendency toward 4 observed in our previous work where the available angular region is limited to | cos θ * | < 0.6 [5]. The fitted value of n = 11.0 ± 0.4 ± 0.4 is much larger than the QCD asymptotic prediction of 6 or 7 [24] but agrees fairly well with n = 10 predicted by a qualitative QCD estimate [55]. For the process γγ → K 0 S K 0 S , according to Refs. [54,55], the W dependence of the cross section is determined by that of the non-leading term in the W region measured by this experiment; the coefficient of the leading term amplitude is much smaller than that of the non-leading term. The results are consistent with the previous analyses [5] with improved statistics and supersede the measurements for the cross section, the χ cJ (1P ) parameters and the slope parameter n.