Amplitude analysis of the $K_{S}K_{S}$ system produced in radiative $J/\psi$ decays

An amplitude analysis of the $K_{S}K_{S}$ system produced in radiative $J/\psi$ decays is performed using the $(1310.6\pm7.0)\times10^{6}$ $J/\psi$ decays collected by the BESIII detector. Two approaches are presented. A mass-dependent analysis is performed by parameterizing the $K_{S}K_{S}$ invariant mass spectrum as a sum of Breit-Wigner line shapes. Additionally, a mass-independent analysis is performed to extract a piecewise function that describes the dynamics of the $K_{S}K_{S}$ system while making minimal assumptions about the properties and number of poles in the amplitude. The dominant amplitudes in the mass-dependent analysis include the $f_{0}(1710)$, $f_{0}(2200)$, and $f_{2}^\prime(1525)$. The mass-independent results, which are made available as input for further studies, are consistent with those of the mass-dependent analysis and are useful for a systematic study of hadronic interactions. The branching fraction of radiative $J/\psi$ decays to $K_{S}K_{S}$ is measured to be $(8.1 \pm 0.4) \times 10^{-4}$, where the uncertainty is systematic and the statistical uncertainty is negligible.

j Currently at: Center for Underground Physics, Institute for Basic Science, Daejeon 34126, Korea (Dated: August 22, 2018) An amplitude analysis of the KSKS system produced in radiative J/ψ decays is performed using the (1310.6±7.0)×10 6 J/ψ decays collected by the BESIII detector. Two approaches are presented. A mass-dependent analysis is performed by parameterizing the KSKS invariant mass spectrum as a sum of Breit-Wigner line shapes. Additionally, a mass-independent analysis is performed to extract a piecewise function that describes the dynamics of the KSKS system while making minimal assumptions about the properties and number of poles in the amplitude. The dominant amplitudes in the mass-dependent analysis include the f0(1710), f0(2200), and f 2 (1525). The mass-independent results, which are made available as input for further studies, are consistent with those of the massdependent analysis and are useful for a systematic study of hadronic interactions. The branching fraction of radiative J/ψ decays to KSKS is measured to be (8.1±0.4)×10 −4 , where the uncertainty is systematic and the statistical uncertainty is negligible.

I. INTRODUCTION
The nature of meson states with scalar quantum numbers has been a topic of great interest for several decades. This is due in part to the expectation that the lightest glueball state should have scalar quantum numbers [1][2][3][4]. Evidence for a glueball state would support longstanding predictions that massive mesons can be generated by gluon self-interactions. Sophisticated studies of experimental data are necessary to observe the effects of gluonic interactions due to the complication of mixing between glueball and conventional quark bound states.
Despite the availability of a large amount of data on ππ and KK scattering in the low mass region, the existence and characteristics of isoscalar scalar (I G J P C = 0 + 0 ++ ) and tensor (0 + 2 ++ ) states remain controversial. The presence of many broad, overlapping states complicates the study of the scalar spectrum, which is poorly described by the most accessible analytical methods [5]. Nonetheless, coupled-channel analyses using the K-matrix formalism have recently produced measurements [6] and dispersive analyses have been directed toward understanding the scalar meson spectrum in the lowest mass region [7]. The BESIII Collaboration has made considerable efforts to improve the knowledge of the scalar and tensor meson sector with a series of amplitude analyses. A mass-dependent (MD) amplitude analysis of radiative J/ψ decays to ηη, using 225 million J/ψ events, describes the scalar spectrum with contributions from the f 0 (1500), f 0 (1710), and f 0 (2100) states [8].
The tensor spectrum appears to be dominated by the f 2 (1525), f 2 (1810), and f 2 (2340) states. BESIII also determined that the f 2 (2340) dominates the tensor spectrum in raditive J/ψ decays to φφ in an amplitude analysis with 1311 million J/ψ events [9]. Additionally, the results of a mass-independent (MI) amplitude analysis of the π 0 π 0 system produced in radiative J/ψ decays include a piecewise function that describes the dynamics of the π 0 π 0 system as a function of invariant mass [10]. These results are useful for developing models that describe hadron dynamics. With the inclusion of additional data from radiative charmonium decays, in particular for the K S K S system, an interpretation of the scalar and tensor meson states may become more clear.
Radiative J/ψ decays to two pseudocalars are a particularly attractive environment in which to study the low mass scalar and tensor meson spectra due to the relative simplicity of an amplitude analysis. Conservation of angular momentum and parity restricts the accessible amplitudes to only those with J P C = even ++ . Radiative J/ψ decays to K + K − have been studied by the MarkIII [11], DM2 [12], and BES [13] Collaborations. The BESII Collaboration performed an amplitude analysis on the K + K − and K S K S system in radiative J/ψ decays, using both a bin-by-bin and global analysis, but the spectrum was limited to less than 2 GeV/c 2 due to the presence of significant backgrounds in the charged channel [14]. A recent comprehensive study of the twopseudoscalar meson spectrum from radiative J/ψ and ψ decays was performed using a 53 pb −1 sample of events at center-of-mass energy √ s = 3.686 GeV taken with CLEO-c [15]. This analysis did not implement a full amplitude analysis, but rather used a Breit-Wigner resonance formalism.
In this paper, we present two independent amplitude analyses of the K S K S system produced in radiative J/ψ decays using the 1311 million J/ψ events collected with the BESIII detector [16]. A MD amplitude analysis parametrizes the K S K S invariant mass spectrum as a coherent sum of Breit-Wigner line shapes, with the goal of extracting the resonance parameters of intermediate states. In addition, a MI amplitude analysis is performed to extract the function that describes the dynamics of the K S K S system using the same method as that described in Ref. [10]. The neutral channel provides a clean environment to study the scalar and tensor meson spectra as it does not suffer from significant backgrounds such as J/ψ → KKπ 0 , which are present in the charged channel J/ψ → K + K − π 0 .

II. BESIII DETECTOR
The BESIII detector is a magnetic spectrometer operating at the Beijing Electron Positron Collider (BEPCII) [17], which is a double ring e + e − collider with center-of-mass energies between 2.0 and 4.6 GeV. The BESIII detector covers a geometrical acceptance of 93% of 4π and consists of a small-celled, helium-based main drift chamber (MDC) which provides momentum and ionization energy loss (dE /dx ) measurements for charged particles; a plastic scintillator time-of-flight system (TOF) which is used to identify charged particles; an electromagnetic calorimeter (EMC), made of CsI(Tl) crystals, that is used to measure the energies of photons and provide trigger signals; and a muon system (MUC) made of resistive plate chambers. A superconducting solenoid magnet provides a uniform magnetic field within the detector. The field strength was 1.0 T during data collection in 2009, but was reduced to 0.9 T during the 2012 running period. The momentum resolution of charged particles is 0.5% at 1.0 GeV/c. The dE /dx measurements provide a resolution better than 6% for electrons from Bhabha scattering. For a 1.0 GeV photon, the energy resolution can reach 2.5% (5%) in the barrel (endcaps) of the EMC, and the position resolution is 6 mm (9 mm). The timing resolution of TOF is 80 ps in the barrel and 110 ps in the endcaps, corresponding to a 2σ K/π separation for momenta up to about 1.0 GeV/c. The spatial resolution of the MUC is better than 2 cm.

III. DATA SETS
This study uses 1311 million J/ψ events collected with the BESIII detector at BEPCII in 2009 and 2012 [16]. An inclusive MC sample of 1225 million J/ψ events generated with the kkmc [18] generator is used for background studies. The main known decay modes are generated using besevtgen [19,20] with branching fractions set to the world average values according to the Particle Data Group (PDG) [5]. The remaining decays are generated according to the Lundcharm model [21].
The K S K S invariant mass distribution of the signal channel in the inclusive MC sample does not resemble that in the data sample. Therefore, for event selection purposes, an exclusive MC sample containing 1 million J/ψ decays to γK S K S is generated according to preliminary results of the MD amplitude analysis. While it does not contain all of the amplitudes in the nominal results of the MD analysis, this MC sample more closely resembles the data and is used to provide a more reliable approximation of the signal to optimize event selection criteria.
An exclusive MC sample, consisting of 5 million J/ψ → γK S K S (K S → π + π − ) events, generated flat in phase space is used for normalization purposes in the MD analysis. A similar exclusive MC sample is used to calculate the normalization integrals for the MI analysis and consists of 110,000 events per 15 MeV/c 2 bin of K S K S invariant mass, with a total of 14.74 million events for the full spectrum. This sample is generated flat in the phase space of each K S K S invariant mass bin, with the result that the overall exclusive MC sample is flat in the distribution of K S K S invariant mass.

IV. EVENT SELECTION CRITERIA
The final state of interest consists of two pairs of charged pions and one photon. Thus the candidate events are required to have at least four good charged tracks whose net charge is zero and at least one good photon. Charged tracks are required to have a polar angle θ that satisfies |cos θ| < 0.93. Each track is assumed to be a pion and no particle identification (PID) restrictions are applied. Each photon is required to have an energy deposited in the EMC greater than 25 MeV in the barrel region (|cos θ| < 0.80) or greater than 50 MeV in the endcaps (0.86 < |cos θ| < 0.92), where θ is the angle between the shower direction and the beam direction, and must fall within the event time (0 ≤ t ≤ 700 ns).
The tracks of each π + π − pair are fitted to a common vertex. Backgrounds that do not contain K S decays are suppressed by restricting L/σ L , where L is the signed flight distance between the common vertex of the π + π − pair and the run-averaged primary vertex, which is taken as the interaction point, and σ L is its uncertainty. For each K S candidate in an event, L/σ L is required to be greater than zero and the value (L1/σ L1 ) 2 + (L2/σ L2 ) 2 is required to be greater than 2.2, where L1 and σ L1 are the distance and uncertainty of one K S , and L2 and σ L2 are those for the other K S in the event.
After the above restrictions are applied, a sixconstraint (6C) kinematic fit is performed to all possible γK S K S combinations, with no charged track used twice in any combination. The 6C kinematic fit consists of four constraints on the energy-momentum of the final state relative to the initial state and one constraint each on the invariant mass of each π + π − pair. The charged track momenta used in the kinematic fit are the updated values after the vertex fit. The χ 2 6C is required to be less than 60. No events have more than one combination of final state particles that survive the above event selection criteria.
A total of 165,137 events survive the event selection criteria. The K S K S and γK S invariant mass spectra are shown in Fig. 1. There are three significant peaks in the K S K S mass spectrum around 1.5, 1.7, and 2.2 GeV/c 2 . The two structures in the γK S spectrum are kinematic reflections from states decaying to K S K S . Figure 2 shows the corresponding Dalitz plots for the data and exclusive MC samples.
The potential backgrounds are studied with the 1225 million J/ψ events of the inclusive MC sample, which is also subjected to the event selection criteria described above. The total amount of backgrounds estimated from the inclusive MC sample is about 0.5% of the size of the data sample. The continuum backgrounds (e + e − → γK S K S without a J/ψ intermediate state) are investigated with a data sample collected at a center-of-mass energy of 3.08 GeV. Only 81 events survive, representing approximately 1,185 events, i.e. 0.7%, of the on-peak data sample after scaling by luminosity and cross section. All backgrounds are ignored in the amplitude analyses.

V. AMPLITUDE ANALYSIS
Amplitude analyses, also called partial wave analyses (PWAs), are typically carried out by modeling the dynamics of particle interactions as a coherent sum of resonances. Such "mass-dependent" (MD) analyses have the benefit that model parameters like Breit-Wigner masses and widths can be related to the properties of the scattering amplitude in the complex s plane, where s is the invariant mass squared of the two-body system. Alternatively, a "mass-independent" (MI) amplitude analysis measures the dynamical amplitude as a function of invariant mass by fitting the sample bin-by-bin while making minimal model assumptions. The results of such an anal- ysis are useful for the development of dynamical models that can be subsequently optimized using experimental data. Each of these methods has benefits and drawbacks as discussed, for example, in Ref. [10]. The correspondence between the model parameters of a MD analysis and the analytic structure of the K S K S amplitude is uncertain due to the presence of broad, overlapping states. On the other hand, the MI analysis suffers from the presence of mathematical ambiguities resulting in multiple sets of optimal parameters in each mass region. The results of the MI analysis are also presented under the assumption of Gaussian errors. This is a necessary step to make the results useful for subsequent analyses, but one that cannot be validated in general. We make use of both analysis methods in this study.
A. MD amplitude analysis

MD amplitude analysis formalism
The MD amplitude analysis is based on the covariant tensor formalism [22]. For radiative J/ψ decays to mesons, the general form of the covariant tensor ampli-tude is where ψ µ (m 1 ) is the J/ψ polarization four-vector, e ν (m 2 ) is the polarization four-vector of the photon and U µν i is the partial wave amplitude with coupling strength determined by a complex parameter Λ i . The U µν i for the intermediate states is constructed from the four-momenta of the daughter particles. The corresponding amplitudes can be found in Ref. [22]. In the MD amplitude analysis, an intermediate resonance is described with the relativistic Breit-Wigner formula with a constant width: where M and Γ are the mass and width of the resonance, respectively, and √ s is the invariant mass of the K S K S system.
Following the convention of Ref. [9], the probability to observe an event characterized by the set of kinematics ξ is where (ξ) is the detection efficiency, ω(ξ) ≡ dσ dΦ is the differential cross section, and dΦ is the standard element of phase space. The full differential cross section is where dξω(ξ) (ξ) ≡ σ is the measured total cross section. A(J P C ) is the full amplitude for all resonances whose spin-parities are J P C . Only K S K S resonances with J P C = 0 ++ , 2 ++ and 4 ++ are considered. For the γK S system, the K 1 and K * resonances are considered. The non-resonant processes are described with a broad resonance whose width is fixed at 500 GeV/c 2 . The complex coupling strength and resonance parameters for each amplitude are determined by an unbinned maximum likelihood fit. The joint probability density for observing N events in the data sample is In practice, the likelihood maximization is achieved by minimizing S = −lnL. The fit is performed based on the GPUPWA framework [23], which takes advantage of parallelization of calculations using Graphical Processing Units (GPUs) to improve computational performance.

MD analysis results
The MD amplitude analysis is performed by assuming the presence of certain expected resonances and then studying the significance of all other accessible resonances listed in the PDG [5]. In Fig. 1, the three structures in the K S K S invariant mass spectrum near 1.5, 1.7, and 2.2 GeV/c 2 indicate the presence of the resonances f 2 (1525), f 0 (1710), and f 0 (2200). These resonances are therefore included in the base solution of the MD analysis. The existence of additional resonances with J P C = 0 ++ , 2 ++ , and 4 ++ above the K S K S threshold and listed in the PDG as well as the intermediate K 1 and K * resonances are then tested. In light of the results of an amplitude analysis of J/ψ decays to φK + K − and φπ + π − by BESII that suggests the presence of an f 0 (1790) [24] that is distinct from the f 0 (1710), the f 0 (1790) is also considered in the MD analysis. The statistical significance of a resonance is evaluated using the difference in log-likelihood, ∆S = − ln L + ln L 0 , and the change in the number of free parameters. Here ln L is the log-likelihood when the amplitude of interest is included and ln L 0 is the log-likelihood without the additional amplitude.
From the set of additional accessible resonances, the one that yields the greatest significance is added to the set of amplitudes if its significance is greater than 5σ. For a wide resonance, the yield must also be larger than 1%. After testing each additional amplitude, the nom- , and f 2 (2340) intermediate states decaying to K S K S as well as the K 1 (1270) and K * (892) intermediate states decaying to γK S . The non-resonant amplitudes for the K S K S system with J P C = 0 ++ and 2 ++ , described by phase space, are also included.
The resonance parameters, i.e. masses and widths, of the dominant 0 ++ and 2 ++ resonances are optimized in the MD analysis. The resonance parameters are listed in Table I, where the parameters listed with uncertainties are optimized while the other parameters are fixed to their PDG values. The systematic uncertainties, which are discussed below, include only those related to the MD analysis. In the resonance parameter optimization, the mass and width of each resonance are optimized by scanning. The values corresponding to the minimum S are taken as the optimized values. The product branching fraction for an intermediate state X is determined according to: where N X is the number of events for the given intermediate state X obtained in the fit, N J/ψ is the total number of J/ψ events, and B K S →π + π − is the branching fraction of K S → π + π − , taken from the PDG [5]. The branching fraction for each process with a specific intermediate state is summarized in Table I. For the decay J/ψ → K S K * (892) with K * (892) → γK S , the measured branching fraction is 6.28 +0.16 −0.17

+0.59
−0.52 × 10 −6 , which is about 3σ away from the product branching fractions taken from the PDG, 10.8 ± 1.2 × 10 −6 . The overall branching fraction for radiative J/ψ decays to K S K S is determined to be (8.29 ± 0.02) × 10 −4 , where the uncertainty is statistical only.
The projections of the K S K S and γK S invariant mass spectra and the angular distributions of the global fit are shown in Fig. 3 and Fig. 4, respectively. The pull distributions of the fit relative to the data are also shown. Given the small statistical uncertainties for such a large data sample, the pulls tend to fluctuate above one. A series of additional checks are also performed for the nominal solution. If the f 0 (1710) and f 0 (1790) are replaced with a single resonance whose mass and width are optimized, S increases by 72.9, indicating that the model of two resonances in this vicinity is preferred over the single resonance model. The f 0 (2200) is also replaced by f 0 (2100) and f 0 (2200) states, but S only decreases by 4.7, corresponding to a significance of less than 5σ. Therefore the parameters for these resonances are set to their PDG values.
In addition to the resonances included in the nominal solution, the existence of extra resonances is also tested. For each additional resonance listed in the PDG, a significance is evaluated with respect to the nominal solution. No additional resonance that yields a significance larger than 5σ also has a signal yield greater than 1% of the size of the data sample. Additionally, an extra f 0 , f 2 , f 4 , K * or K 1 amplitude is included in the fit to test for the presence of an additional unknown resonance. This test is carried out by including an additional resonance in the fit with a specific width (50, 150, 300, or 500 MeV/c 2 ) and a scanned mass in the acceptable region. No evidence for an additional resonance is observed. The scan of the 2 ++ resonance presents a significant contribution around 2.3 GeV/c 2 , with a statistical significance larger than 5σ and a contribution over 1%. However, this hypothetical resonance interferes strongly with the f 2 (2340) due to their similar masses and widths, and is therefore excluded from the optimal solution. B. MI amplitude analysis

MI amplitude analysis formalism
The MI amplitude analysis follows the same general procedure as that described in Ref. [10]. The amplitudes are extracted independently in bins of K S K S invariant mass. Only the 0 ++ and 2 ++ amplitudes are found to be significant in the analysis. Under the inclusion of a 4 ++ amplitude, no bins yield a difference in S equivalent to a 5σ difference. Only one bin yields such a difference for the case of a K * K 0 amplitude, where the K * decays to γK 0 . The K * K 0 amplitude is spread over many K S K S bins and therefore does not contribute significantly to any individual K S K S invariant mass bin. The effect on the results for the case of a possible additional amplitude is taken as a systematic uncertainty. The amplitudes for radiative J/ψ decays to K S K S are identical to those for radiative J/ψ decays to π 0 π 0 . In brief, the amplitude is constructed as where x = {θ γ , φ γ , θ K , φ K } is the position in phase space, s is the invariant mass squared of the K S K S pair, M is the polarization of the J/ψ, and λ γ is the helicity of the radiative photon. Here, both M and λ γ may have values of ±1. The amplitude is then factorized, with one piece describing the radiative transition to an intermediate state X and the other describing the strong interaction dynamics  [5]. The branching fractions and significance for each resonance is also given. When two uncertainties are given for a branching fraction, the first and second uncertainties are statistical and systematic, respectively. The systematic uncertainties due to overall normalization affect the branching fractions, but have little effect on the mass and width parameters. where j is the angular momentum of the intermediate state and J γ indexes the radiative multipole transitions. Any pseudoscalar-pseudoscalar final states that may rescatter into the K S K S final state are accounted for in the sum over X. In the radiative multipole basis, the amplitudes include an E1 component for J P C = 0 ++ and E1, M 2, and E3 components for J P C = 2 ++ . Finally, the amplitude may be written where V j,Jγ (s) is the coupling to the state with characteristics j and J γ . This coupling factor includes the complex function that describes the K S K S dynamics as well as the coupling for the radiative decay, which cannot be separated. The piece of the amplitude that describes the angular distributions, A M,λγ j,Jγ ( x), is determined by the kinematics of an event.
In the MI analysis, the data sample is binned as a function of K S K S invariant mass, under the assumption that the part of the amplitude that describes the strong interaction dynamics is constant over a small range of s, This is done to avoid making strong model dependent assumptions about the dynamical function. The couplings are then taken as free parameters in an extended maximum likelihood fit in each mass bin. In this way, a table of complex numbers is extracted representing the free parameters in each bin that describe the K S K S interaction dynamics.
The density of events at some position in phase space x is given by the intensity function, where the free parameters are constrained to be the same for each piece of the incoherent sum over the (unmeasured) observables of the interaction. The observables include the polarization of the J/ψ, M = ±1, and the helicity of the radiative photon, λ γ = ±1.
The intensity for the amplitude in bin k, bounded by s k and s k+1 , indexed by j and J γ is given by where the fit parameters, V k j,Jγ , are the product of V k j,Jγ and the square root of the size of the phase space in bin k. The intensities presented in Figs. 5 and 6 as well as in the supplemental materials [29] for the MI analysis are corrected for detector acceptance and efficiency.

Ambiguities
The MI amplitude analysis is complicated by the presence of ambiguities. A phase convention is applied to remove trivial ambiguities created by the freedom to rotate the overall amplitude by π or to reflect it over the real axis in the complex plane. This freedom comes from the fact that the intensity is constructed from a sum of absolute squares. Non-trivial ambiguities are discussed in detail in Ref. [10] and are due to the possibility for amplitudes with the same quantum numbers to have different phases. As shown in Ref. [10], only two ambiguous solutions are present for the case of J/ψ radiative decays to two pseudoscalars if only the 0 ++ and 2 ++ amplitudes are considered. Both solutions are presented for bins in which the ambiguous solutions are not degenerate. If additional amplitudes are introduced, the number of am-biguities would increase.

MI analysis results
The intensities for each amplitude and the phase differences relative to the reference amplitude, 2 ++ E1, are plotted in Fig. 5 and Fig. 6, respectively. Several bins exhibit two ambiguous solutions, but for many bins, the ambiguous partner is degenerate. An arbitrary phase convention is applied in which the phase difference between the 0 ++ and 2 ++ E1 amplitudes is required to be positive. For much of the spectrum, the ambiguous solutions do not exhibit two distinct continuous sets of solutions though there is some indication that two distinct sets of solutions exist below about 1.5 GeV/c 2 .
Finally, the branching fraction for radiative J/ψ decays to K S K S is determined according to Here, N γK S K S is the acceptance corrected signal yield determined by summing the total intensity from each K S K S invariant mass bin in the MI analysis results, N bkg is the acceptance corrected background contamination determined from the inclusive MC and continuum data samples, and N J/ψ is the total number of J/ψ events in the data sample. An efficiency correction γ is applied in order to extrapolate the K S K S invariant mass spectrum down to a radiative photon energy of zero and is determined by calculating the fraction of phase space that is removed by restricting the energy of the radiative photon. This extrapolation results in an increase in the total number of events by 0.02%, so γ is taken to be 0.9998.
To determine N bkg , the efficiency correction for the inclusive MC background and continuum samples is assumed to be the same as that for the data sample. That is, N bkg is determined according to where N γK S K S ,k is the acceptance corrected signal yield in bin k, N acc γK S K S ,k is the number of events in the data sample for bin k, and N mc,k is the number of background events in bin k according to the inclusive MC and continuum samples. This method gives a background fraction, N bkg /N γK S K S , of about 1.14%, which is roughly consistent with the approximation of a background contamination of 1.11% according to the number of background events in the inclusive MC sample relative to the size of the data sample. According to Eq. (13), the branching fraction for radiative J/ψ decays to K S K S is determined to be (8.10 ± 0.02) × 10 −4 , where only the statistical uncertainty is given.
It is also important to note that the MI analysis results are only valid in the Gaussian limit. As discussed in the amplitude analysis of J/ψ decays to γπ 0 π 0 [10], this assumption cannot be guaranteed for all parameters in the analysis. Therefore, the use of these results may not produce statistically rigorous values for parameters of interest. Rigorous values of model parameters can only be reliably extracted by fitting a model directly to the data.

C. Discussion
The nominal results of the MI and MD analyses are in good agreement. A comparison of the total 0 ++ and 2 ++ intensities without acceptance correction are shown in Fig. 7. The results of the MI analysis show significant features in the 0 ++ amplitude just above 1.7 GeV/c 2 and just below 2.2 GeV/c 2 , consistent with the f 0 (1710) and f 0 (2200), respectively. The former of these states is often cited as a scalar glueball candidate [25,26]. Additional structure above 2.3 GeV/c 2 suggests the need for another state in this region. This is in agreement with the MD analysis, which suggests that the f 0 (1710) and f 0 (2200) dominate the scalar spectrum and also includes an f 0 (2330). Additionally, the scalar spectrum near and below 1.5 GeV/c 2 shows a complicated structure. The presence of the f 0 (1370) and f 0 (1500) may be necessary to describe this region, as in the MD results.
The 2 ++ amplitude extracted in the MI analysis is dominated by a structure near 1.5 GeV/c 2 , which may reasonably be interpreted as the f 2 (1525), in agreement with the MD analysis and Ref. [14]. The 2 ++ amplitude near 1.2 GeV/c 2 in the MI results suggests the presence of a state like the f 2 (1270) as in the MD analysis.
The branching fraction for the MD analysis does not take into account the small remaining backgrounds. Therefore, the branching fraction measurement from the MI analysis is taken as the nominal result. The measurement is also repeated for the MI analysis without subtracting the backgrounds. The result is (8.20 ± 0.02) × 10 −4 . The difference between this value and that determined in the MD analysis is taken as a systematic uncertainty. The small discrepancy is likely due to the difference in the efficiency calculation for the two methods. The efficiency for the MD analysis depends on the fitting result so the fit quality can have an influence on the branching fraction. The branching fraction measurement is dominated by systematic effects, which are discussed below.

VI. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties for this analysis are divided into three different categories. The first is the systematic uncertainty due to the overall normalization of the results. Sources of this type of uncertainty include the K S reconstruction, the 6C kinematic fit, and the photon detection efficiency, which are described in Section VI A. Additional sources of uncertainty related to the overall normalization include the total number of J/ψ events, which is taken from Ref. [16], the decay branching fraction of K S → π + π − , the analysis method, and the remaining backgrounds. The systematic uncertainties related to the overall normalization are described in detail in Section VI A and summarized in Table II. The other sources of systematic uncertainty are specific to the MD or MI analysis methods and are described in Sections VI B and VI C.

A. Systematic uncertainty related to the overall normalization
The K S reconstruction efficiency is studied with a control sample of J/ψ → K * ± (892)K ∓ events, where K * ± (892) → K S π ± . A fit is applied to the missing mass squared recoiling against the K ± π ∓ system to determine the fraction of candidate events that pass the K S selection requirements given above. In the fit, the signal shape is taken from an exclusive MC sample, convolved with a Gaussian function. The background is fixed to the shape of the backgrounds extracted from the inclu-sive MC sample. The momentum weighted difference in the K S reconstruction efficiency between the data and MC samples is taken as the associated systematic uncertainty. The total uncertainty due to K S reconstruction for the event topology of interest is determined to be 4.1%.
A control sample of ψ → γχ c0,2 , with χ c0,2 → K S K S is used to estimate the uncertainty associated with the 6C kinematic fit. The efficiency is the ratio of the signal yields with and without the kinematic fit requirement, χ 2 6C < 60. The difference in efficiency between the data and MC samples, 1.2%, is taken as the systematic uncertainty.
The photon detection efficiency of the BESIII detector is studied using a control sample of J/ψ decays to π + π − π 0 , where the π 0 decays into two photons. The largest difference in the photon detection efficiency for the inclusive MC sample with respect to that for the data sample is taken as the systematic uncertainty due to photon reconstruction. The systematic uncertainty is determined to be 0.5% for photons with an angular distribution of |cos θ| < 0.8 and 1.5% for photons that fall in the endcap region (0.86 < |cos θ| < 0.92). For radiative J/ψ decays to K S K S , 93% of the reconstructed photons fall in the barrel region. Therefore, the systematic uncertainty due to the photon detection efficiency is determined to be 0.6%.
The amplitude analyses are performed under the assumption of no backgrounds. Therefore, an uncertainty due to the background events is assigned. Conservative systematic uncertainties equal to 100% of the background contamination are attributed to each of the inclusive MC and continuum background types. The systematic uncertainty associated with the remaining backgrounds is about 0.5% for the backgrounds from the inclusive MC sample and about 0.7% for the continuum backgrounds.
The difference in the branching fraction for radiative J/ψ decays to K S K S between the MD and MI analyses is taken as a systematic uncertainty due to the analysis method. Both methods are used to determine the branching fraction in the case where background contamination is ignored, yielding a difference of 1.1%.
The total systematic uncertainty for the overall normalization is determined by assuming all of the sources described above are independent. The individual uncertainties are combined in quadrature, resulting in an uncertainty of 4.6%.

B. Systematic uncertainties related to the MD analysis
Uncertainties due to possible additional amplitudes in the MD analysis are estimated by adding, individually, the most significant amplitudes from the extra resonance checks described above. These additional amplitudes include the K 1 (1400), K * PHSP, f 0 (2100), f 2 (1810) and 4 ++ PHSP. The changes in the measurements relative to  the nominal results are taken as systematic uncertainties.
In the optimal solution of the MD analysis, the resonance parameters of some amplitudes are fixed to PDG values [5]. An alternative fit is performed in which those resonance parameters are varied within one standard deviation. The changes in the measurements are taken as systematic uncertainties.
In addition to the global uncertainty due to the K S reconstruction efficiency, an uncertainty related to the difference in the momentum dependence of the K S reconstruction efficiency between the data and MC simulation is considered in the MD analysis. The reconstruction efficiency of the phase space MC sample used in the MD analysis is corrected and the fit is repeated with the nominal central values. The differences in the branching fraction measurements between these and the nominal results are taken as systematic uncertainties.
For some parameters, the systematic variations leave the central value unchanged, indicating that the systematic uncertainty is negligible. The total systematic uncertainties related to the MD analysis are given in Table I.

C. Systematic uncertainties related to the MI analysis
The Dalitz plot in Fig. 1 (a) shows a K * K 0 amplitude, where the K * decays to γK 0 , especially for the high K S K S mass region. This amplitude is also apparent in the MD analysis results. In the MI analysis, the K * K 0 amplitude is spread over many K S K S invariant mass bins and does not contribute significantly in any individual mass bins. With the inclusion of a K * K 0 amplitude, the results of the MI analysis do not change significantly. This suggests that the MI analysis is not sensitive to the K * K 0 amplitude, so it is neglected.
Only amplitudes with J P C = even ++ are allowed in radiative J/ψ decays to K S K S . The results of the MD analysis and the nominal results of the MI analysis only include 0 ++ and 2 ++ amplitudes and no 4 ++ amplitude. Under the inclusion of a 4 ++ amplitude, the results of the MI analysis do not change significantly. This suggests that the 4 ++ amplitude does not contribute or that the MI analysis is not sensitive to it, if it does exist.
A study of the effect that an additional 4 ++ amplitude would have on the MI analysis suggests that deviations occur on the order of the statistical uncertainties of the data sample [10]. Therefore, the systematic uncertainty due to the effect of ignoring a possible additional amplitude is estimated to be of the same order as the statistical uncertainties of the MI results.

VII. CONCLUSIONS
An amplitude analysis of the K S K S system produced in radiative J/ψ decays has been performed using two complementary methods. A mass-dependent amplitude analysis is used to study the existence and coupling of various intermediate states including light isoscalar resonances. The dominant scalar amplitudes come from the f 0 (1710) and f 0 (2200), which have production rates in radiative J/ψ decays consistent with predictions from lattice QCD for a 0 + glueball and its first excitation [27]. The production rate of the f 0 (1710) is about one order of magnitude larger than that of the f 0 (1500), which suggests that the f 0 (1710) has a larger overlap with the glueball state compared to the f 0 (1500). The tensor spectrum is dominated by the f 2 (1525) and f 2 (2340). Recent Lattice QCD predictions for the production rate of the pure gauge tensor glueball in radiative J/ψ decays [28] are consistent with the large production rate of the f 2 (2340) in the K S K S , ηη [8], and φφ [9] spectra.
The mass-dependent results are consistent with the results of a mass-independent amplitude analysis of the K S K S invariant mass spectrum. The mass-independent results are useful for a systematic study of hadronic interactions. The intensities and phase differences for the amplitudes in the mass-independent analysis are given in supplemental materials [29]. A more comprehensive study of the light scalar meson spectrum should benefit from the inclusion of these results with those of similar reactions. Details concerning the use of these results are given in Appendix C of Ref. [10].
Finally, the branching fraction for radiative J/ψ decays to K S K S is determined to be (8.1±0.4)×10 −4 , where the uncertainty is systematic and the statistical uncertainty is negligible.