text:
# Statistical Analysis of Different 𝜈̄_𝜇→𝜈̄_𝑒 Searches

## I Introduction

Currently, there are three neutrino anomalies interpreted as evidence for neutrino oscillations, namely the atmospheric atmo , the solar solar  and the LSND LSNDfinal  anomaly with three distinct squared neutrino mass differences. However, oscillations between the three Standard Model (SM) neutrinos are described by only two independent neutrino mass–squared differences, allowing only two of the above anomalies as being due to oscillations in a SM scenario extended by massive neutrinos. Before enlarging the neutrino sector by some additional, unknown sterile neutrinos, each of the experimental results has to be checked independently.

Over the last years, the controversial results of the two experiments LSND (Liquid Scintillator Neutrino Detector at LANSCE, Los Alamos, USA) and KARMEN (KArlsruhe Rutherford Medium Energy Neutrino experiment at ISIS, Rutherford, UK) both searching for neutrino oscillations $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ in a short baseline regime, have led to intense discussions. The two experiments are similar as they use $`\overline{\nu }_\mu `$ beams from the $`\pi ^+`$-$`\mu ^+`$ decay at rest (DAR) chain $`\pi ^+`$$``$ $`\mu ^+`$\+ $`\nu _\mu `$ followed by $`\mu ^+`$ $``$ e+ \+ $`\nu _e`$ \+ $`\overline{\nu }_\mu `$ with energies up to 52 MeV. Furthermore, both experiments are looking for $`\overline{\nu }_e`$ from $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations according to the two–flavor oscillation probability 333In a comparison of different oscillation searches, one generally has to use a complete 3– or 4–dimensional description of neutrino masses and mixings, leading to the general oscillation probability. Since both experiments under investigation here have similar parameters and, most important, search for oscillations in the same appearance mode, the simplified formula holds for the purpose of the following analysis.

$$P(\overline{\nu }_\mu \overline{\nu }_e)=\mathrm{sin}^2(2\mathrm{\Theta })\mathrm{sin}^2\left(1.27\frac{\mathrm{\Delta }m^2L}{E_\nu }\right)$$
(1)
with the difference of the squared masses $`\mathrm{\Delta }m^2`$ in eV2/c4, the flight length of the neutrino $`L`$ in meters and the neutrino energy $`E_\nu `$ in units of MeV. The detection reaction p ( $`\overline{\nu }_e`$ , e+ ) n  provides a spatially correlated delayed coincidence signature of a prompt e+ and a subsequent neutron capture signal.

This paper describes a combined statistical analysis of the final LSND and KARMEN 2 results. There are significant changes to the treatment of LSND data in the final LSND paper LSNDfinal . The present paper uses only the decay-at-rest (DAR) $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ data, and leaves out the impact of the decay-in-flight (DIF) $`\nu _\mu `$$``$$`\nu _e`$ data. As shown later, the LSND data has a favored region in the 7 eV2/c4 region of $`\mathrm{\Delta }m^2`$ as a consequence of this change. This is in contrast to the LSND result LSNDfinal  which disfavors this solution at 90% C.L. because of the inclusion of the DIF flux. A motivation for this change is that a direct comparison of the two experiments is less model dependent if one considers the LSND DAR flux only, since KARMEN does not have a DIF neutrino flux present in its source.

The analysis is based on a frequentist approach following the suggestions of feldcous . For both experiments, the data are analysed with a maximum likelihood analysis followed by the extraction of confidence levels in a unfied approach. It is not the task nor the purpose of this analysis to treat the apparent disagreement between the two experiments. Rather, assuming no serious systematical error in either experiment or their interpretation with respect to $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillation, we apply a consistent statistical analysis to quantitatively establish their level of compatibility. The work follows an earlier analysis NJP  based on intermediate data sets.

There are two separate questions investigated in this paper that relate to the LSND and KARMEN data. The first question is: at what level are the two data sets compatible with each other? This question addresses whether or not the two experiments force one to draw opposite conclusions from their respective data.

If one assumes that the data sets are compatible, the second question becomes: what are the most likely regions for the oscillation parameters? This method is somewhat analogous to treating the KARMEN data as coming from a ’near detector’ at 17 meters, and treating the LSND data as coming from a ’far detector’ at 30 meters. It ignores the question of systematic differences between the experiments, such as whether or not, for example, an unaccounted source of $`\overline{\nu }_e`$ background neutrinos was present in the LSND data and not present in the KARMEN data.

The paper is organized as follows: After a short overview of the experimental setups of both experiments, section II describes the data analyses leading to likelihood functions of the oscillation parameters. Section III is devoted to the extraction of confidence regions for the individual experiments. In section IV, the experimental results are combined statistically, extracting levels of compatibility as well as oscillation parameters compatible with both experiments. We conclude with implications of this joint analysis, a comparison with other oscillation searches and an outlook to upcoming experimental tests of the favored oscillation parameters presented in this analysis.

### I.1 The LSND experiment

The source of neutrinos for the LSND experiment was the interaction of the 798 MeV proton beam at the Los Alamos Neutron Science Center (LANSCE), in which a large number of pions, mostly $`\pi ^+`$ are produced. The $`\pi ^{}`$ are mainly absorbed and only a small fraction decay to $`\mu ^{}`$, which in turn are largely captured. Thus, the resulting neutrino source is dominantly due to $`\pi ^+`$$``$ $`\mu ^+`$\+ $`\nu _\mu `$ and $`\mu ^+`$ $``$ e+ \+ $`\nu _e`$ \+ $`\overline{\nu }_\mu `$ decays, most of which decay at rest (DAR). With a very small contamination of $`\overline{\nu }_e/\overline{\nu }_\mu 810^4`$, a measurement of the reaction $`\overline{\nu }_e`$ \+ p $``$ n + $`e^+`$  provides a sensitive way to search for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations. Such events are identified by detection of both the $`e^+`$ and the 2.2 MeV $`\gamma `$ from the reaction p ( n,$`\gamma `$ ) d. In addition, the $`\nu _e`$ flux from $`\pi ^+`$ and $`\mu ^+`$ decay-in-flight (DIF) is very small, which allows a search for $`\nu _\mu `$$``$$`\nu _e`$ oscillations via the measurement of electrons above the Michel electron endpoint from the reaction 12C ( $`\nu _e`$ , e- ) 12N.

The LSND experiment took data over six calendar years (1993-1998). During this period the LANSCE accelerator operated for 17 months, delivering 28 896 C of protons on the production target. The LSND detector bigpaper  consisted of an approximately cylindrical tank 8.3 m long by 5.7 m in diameter. The center of the detector was located 30 m from the beam stop neutrino source. The tank was filled with liquid scintillator consisting of mineral oil and 0.031 g/l of b-PBD. This low scintillator concentration allowed the detection of both Čerenkov light and scintillation light. Photomultiplier time and pulse-height signals were used to reconstruct the track with an average RMS position resolution of $`14`$ cm and an energy resolution of $`7`$ % at the Michel endpoint of 52.8 MeV.

### I.2 The KARMEN experiment

The KARMEN experiment used as neutrino source the pulsed spallation neutron source ISIS of the Rutherford Appleton Laboratory delivering 800 MeV protons. In contrast to the LANSCE source, the protons are extracted from the synchrotron as an intense but narrow double pulse, consisting of two parabolic pulses with a width of 100 ns separated by a peak-to-peak gap of 325 ns. The unique time structure of the ISIS proton pulses allowed a clear separation of neutrino-induced events from any beam unrelated background. The intrinsic contamination of the ISIS $`\overline{\nu }_\mu `$–beam is $`\overline{\nu }_e/\overline{\nu }_\mu =6.410^4`$ Bur96 .

The KARMEN detector drexlin  was a rectangular high resolution liquid scintillation calorimeter, located at a mean distance of 17.7 m from the ISIS target at an angle of 100 degrees relative to the proton beam. The liquid scintillator was enclosed by a multilayer active veto system and a 7000 t steel shielding. The 65 m3 of liquid scintillator consisted of a mixture of paraffin oil ($`75\%`$vol.), Pseudocumene ($`25\%`$vol.) and 2 g/l scintillating additive 1-phenyl-3-mesityl-2-pyrazoline. The liquid scintillator volume was optically separated into 512 independent central modules. Gadolinium was implemented between the module walls for an efficient detection of thermal neutrons Gd ( n,$`\gamma `$ ) with on average 3 $`\gamma `$’s of energy $`E_\gamma =8`$ MeV. The KARMEN detector as liquid scintillator calorimeter was optimized for high energy resolution of $`\sigma _E=11.5\%/\sqrt{E(\text{MeV})}`$. Each event had energy, time and position information, as well as the number of addressed modules and their relative time differences.

The KARMEN 2 experiment took data from February 1997 to March 2001. During this time, protons equivalent to a total charge of 9425 Coulombs have been accumulated on the ISIS target. This corresponds to $`2.7110^{21}`$ $`\overline{\nu }_\mu `$ from the ISIS beam stop target.

## II Data analysis and results as likelihood functions

For both experiments, a $`\overline{\nu }_e`$ signal from $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations consists of a spatially correlated delayed (e+,n) sequence from p ( $`\overline{\nu }_e`$ , e+ ) n . Due to the different experimental techniques, the parameters to identify such sequences vary. Details of the LSND event reconstruction can be found in LSNDdet  and bigpaper , the requirements for event sequences in KARMEN are described in K2paper .

In LSND, a (e+,n) sequence requires a prompt ’electron-like’ event with energy $`E_e>20`$ MeV followed by a low energy $`\gamma `$–event. The information about the delayed event is characterized via a likelihood ratio $`R_\gamma `$: If within one millisecond after the initial event another event is recorded at distance $`\mathrm{\Delta }r250`$ cm, the likelihood ratio $`R_\gamma `$ in energy (PMT hits), time and distance of being a correlated p ( n,$`\gamma `$ ) over an accidental coincidence is calculated, otherwise $`R_\gamma =0`$. Requiring a high likelihood ratio $`R_\gamma `$ selects (e+,n) correlated events with low uncorrelated background, the so-called “gold-plated” events. However, to determine the oscillation parameters $`\mathrm{sin}^2(2\mathrm{\Theta })`$ and $`\mathrm{\Delta }m^2`$ in a likelihood analysis, an event sample with a much looser cut in $`R_\gamma `$ is used, leading to higher efficiency for the oscillation channel, but naturally increasing the background. For the LSND sample analysed here, all the cuts described in LSNDfinal  have been applied with two exceptions: The energy of the prompt event has been restricted to $`20<E_e<60`$ MeV instead of $`20<E_e<200`$ MeV and we applied a loose cut of $`R_\gamma >10^5`$ instead of no $`R_\gamma `$-cut at all. The reason for these changes is twofold. First, the upper energy cut $`E_e<60`$ MeV excludes any significant influence of the $`\nu _\mu `$$``$$`\nu _e`$ channel on the much more sensitive $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ channel. For any value of the oscillation parameter $`\mathrm{\Delta }m^2`$, the expected contribution from $`\nu _\mu `$$``$$`\nu _e`$ remains below 2 % of the $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ signal. The second reason is a more technical one. Applying a unified frequentist approach to the likelihood function gained from the reduced data implies simulation and analysis of a large number of event samples analog to the experiment. Therefore, one is interested not to have too large samples containing mainly background events. With the cuts above, the original event sample is reduced by a factor of 5.5 to 1032 candidates, with an efficiency for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ at high $`\mathrm{\Delta }m^2`$ reduced to 64.4 % of the original one (only 3.2 % of the expected $`\nu _\mu `$$``$$`\nu _e`$ events remain).

Figure 1 shows the event sample comprising 1032 beam-on events. Four variables are used to categorize the events:

The energy of the primary electron $`E_e`$, its spatial distribution along the detector axis $`z`$ 444The spatial coordinate $`z`$ along the detector axis can be easily referred to the neutrino flight path $`L`$ by adding the distance source–center of detector tank. and the angle $`\mathrm{cos}\theta _\nu `$ between the direction of the incident neutrino and the reconstructed electron path. The fourth variable is the likelihood ratio $`R_\gamma `$ for a (e+,n) coincidence. The binning of this distribution with the upper edge of bin number $`N`$ representing a value $`R_\gamma =10^{N/31}`$ has been chosen for practical reasons. Note that the spectra in Figure 1 are projections of a 4-dim space of correlated parameters for each event. Superimposed to the data are shown the beam-unrelated and beam-related background contributions amounting to a total number of $`N_{bgd}^{tot}=971.9`$ events. The expected background contributions are broken down in Table 1, for further details we refer to LSNDfinal .

Formally, each beam-on event $`j`$ of the $`N_{beamon}=1032`$ candidates is assigned a probability $`p_j(\stackrel{}{x})`$ equal to a sum of probabilities $`q_i(\stackrel{}{x})`$ from the backgrounds plus oscillations. The vector $`\stackrel{}{x}`$ describes herein the spectral parameters of each event $`j`$, $`\stackrel{}{x}_j=(E_{ej},R_{\gamma j},\mathrm{cos}\theta _{\nu j},z_j)`$. It then remains to add the $`q_i`$ with expected fractional contributions $`r_i`$ and take the product over all the beam-on events. The likelihood is thus

$$=(\underset{j=1}{\overset{N_{beamon}}{}}p_j),$$
(2)
where

$$p_j(\stackrel{}{x}_j)=\underset{i=1}{\overset{N_{contr.}}{}}q{}_{i}{}^{}(\stackrel{}{x}_j)r_i.$$
(3)
Additionally, two normalization requirements must hold:

$$\underset{i=1}{\overset{N_{contr.}}{}}r_i=1,$$
(4)
and

$$𝑑E_e𝑑R_\gamma d(\mathrm{cos}\theta _\nu )𝑑zq_i(E_e,R_\gamma ,\mathrm{cos}\theta _\nu ,z)=1$$
(5)
for each contribution, $`i`$. Together, these requirements ensure that every observed beam-on event has a probability of occurrence equal to 1.

To allow for the fact that the backgrounds are known to some limited accuracy only, the background is varied by calculating the above likelihood at each point in the ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) plane many times, varying over the expected $`\sigma `$ for each background. For each background configuration, the $``$ is weighted with a Gaussian factor for each background that is off its central value. The background variations are performed in a simplified manner, i.e. the beam-unrelated background (BUB) varies independently and all the beam-related backgrounds (BRBs) are locked together. Finally, the likelihood can be expressed as

$$=𝒟N_{bgd}\mathrm{exp}((N_{bgd}N_{bgd,exp})^2/2\sigma ^2)\underset{j=1}{\overset{N_{beamon}}{}}p_j,$$
(6)
where the $`𝒟N_{bgd}`$ represents, schematically, the background variation described above. Integrating over different background contributions eventually reduces the free parameters of the likelihood procedure from $`N_{contr.}+1`$ to the 2 free parameters of real interest within this analysis, the oscillation parameters:

$$(r_1,\mathrm{}r_{N_{contr.}1},\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)$$
(7)
Figure 2 shows the logarithm of the event–based likelihood

function, $`\mathrm{ln}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)`$, for the 1032 events investigated. The maximum is reached at a parameter combination ($`\mathrm{sin}^2(2\mathrm{\Theta })`$=0.85, $`\mathrm{\Delta }m^2`$=0.055 eV2/c4) corresponding to a total oscillation signal of $`N_{osc}=N(\overline{\nu }_\mu \overline{\nu }_e)+N(\nu _\mu \nu _e)=66.4+0.3=66.7`$. The difference to the no–oscillation hypothesis as expressed in logarithmic likelihood units is

$$\mathrm{ln}(0.85,0.055\text{eV}\text{2}\text{/c}\text{4})\mathrm{ln}(0,0)=23.5$$
(8)
underlining the significance of the additional $`\overline{\nu }_e`$ signal among the event sample. Note that the extracted oscillation events $`N_{osc}`$ scale with the efficiency of the applied cuts and are in good agreement with the best fit result stated in LSNDfinal .

KARMEN 2 collected data from February 1997 through March 2001 corresponding to 9425 C accumulated proton charge on the ISIS target. A spatial coincidence between the initial e+ and the neutron capture of 1.3 m3 was required. Applying all cuts to the data K2paper , 15 (e+,n) candidate sequences were finally reduced. Figure 3 shows the remaining sequences in the appropriate energy, time and spatial windows.

The background components are also given with their distributions. All components except the intrinsic $`\overline{\nu }_e`$ contamination are measured online in different time and energy windows (see Table 2).

The extracted number of sequences is in excellent agreement with the background expectation, consistent with no oscillation signal. To also include the detailed spectral information of each individual event, a maximum likelihood method is applied.

Analogously to the LSND likelihood function, we can define the combined likelihood for the KARMEN sample as

$$=(\underset{j=1}{\overset{N_{beamon}}{}}p_j),$$
(9)
where

$$p_j(\stackrel{}{x}_j)=\underset{i=1}{\overset{N_{contr.}}{}}q{}_{i}{}^{}(\stackrel{}{x}_j)r_i$$
(10)
with $`N_{beamon}=15`$ and $`N_{contr.}=5`$, the oscillation signal as well as 4 background contributions. For KARMEN 2, the vector $`\stackrel{}{x}`$ has, of course, a different definition, namely $`\stackrel{}{x}=(E_{pr},T_{pr},E_{del},\mathrm{\Delta }T,\mathrm{\Delta }x)`$ describing the energy and time relative to beam-on-target of the prompt event, the energy of the delayed event as well as the time difference and spatial correlation of the prompt and delayed event. Again, normalizations (4) and (5) with the appropriate event parameters have to hold. Differing from the likelihood definition in the LSND analysis, for KARMEN with its very small background, all background components are combined into one contribution, therefore reducing the above sum (10) to only the oscillation signal and the total background, $`N_{contr.}=2`$. Furthermore, the weight for the background variation, defined in (6) as a Gaussian factor, is taken as Gamma function, reflecting a continuous variation in a Poisson statistics K2paper .

Figure 4 shows the logarithm of the likelihood as a function of the oscillation parameters ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$). The maximum is reached near the physical boundary $`\mathrm{sin}^2(2\mathrm{\Theta })=0`$, i.e. the no–oscillation case. The sharp drop in $`\mathrm{ln}`$ towards larger values of $`\mathrm{sin}^2(2\mathrm{\Theta })`$ demonstrates that there is no indication for an oscillation signal.

A typical approach to deduce confidence regions from likelihood functions would consist of defining the contour of an area of given confidence by cutting the likelihood function at the corresponding value below the maximum. For a Gaussian shaped likelihood function of two uncorrelated free parameters, the corresponding values would be

$$\mathrm{\Delta }\mathrm{ln}_{\mathrm{Gauss}}=\mathrm{ln}(1\delta /100)=2.3(4.6)$$
(11)
for $`\delta =90(99)\%`$ confidence, respectively Lyons . However, this method leads to correct coverage only in the Gaussian case. As can be seen from Figures 2 and 4 the likelihood functions have distinct features: An oscillatory behavior as a function of the free parameters, the maximum at the physical boundary (for KARMEN), or numerous side maxima with comparable likelihood values (for LSND). In addition, the parameter space in ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) is, in principle, not limited, its metric in terms of prior probability density not unequivocally defined.

## III Individual confidence regions

The method we apply in the following is based on a unified frequentist approach suggested by feldcous  which eliminates the bias that occurs when one decides, after analysing the data, between using a confidence interval (having a positive signal) or an upper confidence limit (having a result compatible with the background). When first presented, it was argued that the suggested ordering principle near the physical boundary should be modified Giunti ,Byron , or that a Bayesian approach to the extracted likelihood function would be more appropriate Dagos . However, many experimental results have been analysed since using the unified approach CWCL , also being described as a standard procedure by Pdg2000 .

The basic idea of the unified frequentist approach is to create a large number of event samples in perfect analogy to an experiment. These samples are created by Monte Carlo using the full event information for the likelihood procedure. For an oscillation hypothesis $`H`$ with given parameters ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$)H, a large sample of Monte Carlo simulations of so-called toy experiments is created. These simulations are based on the detailed knowledge of all experimental resolution functions and the spectral information on the individual background contributions. In addition, they comprise the expected experimental signal for the given oscillation hypothesis 555For each individual sample, the number of background events as well as the number of signal events vary according to the error of the expectation values.. For each likelihood function of the analysed samples, the estimator $`\mathrm{\Delta }\mathrm{ln}=\mathrm{ln}_{max}\mathrm{ln}_H`$ is calculated by comparing each sample’s likelihood maximum with the value for the given input parameters ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$)H. The hypothesis $`H`$ is then accepted at a confidence $`\delta `$, if the estimator $`\mathrm{\Delta }\mathrm{ln}`$ of the experimental likelihood function is contained within the smallest $`\delta `$ % of the simulated estimator distribution. Since, in principle, the estimator distribution itself is a function of the oscillation parameters, a complete statistical analysis consists of a scan of the entire parameter space ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) to extract the according region of confidence.

### III.1 KARMEN

To deduce confidence regions from the KARMEN 2 likelihood function, the parameter space ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) has been scanned extracting the estimator distribution for 8000 different oscillation hypotheses, including the no–oscillation scenario. From these distributions, each containing the information of about 4000 simulated and analysed samples, regions of different confidence levels can be calculated.

Figure 5 shows the regions of the oscillation parameters ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) for confidence levels of 68 %, 90 % and 95 %. Parameter combinations to the right contain the complementary confidence. The contours of the confidence areas are often referred to as exclusion curves of a given confidence. For values $`\mathrm{\Delta }m^2100`$eV2/c4, the contours cut at oscillation amplitudes of $`\mathrm{sin}^2(2\mathrm{\Theta })=0.910^3`$ (68 % C.L.), $`\mathrm{sin}^2(2\mathrm{\Theta })=1.710^3`$ (90% C.L.) and $`\mathrm{sin}^2(2\mathrm{\Theta })=2.210^3`$ (95% C.L.).

The KARMEN 2 sensitivity, defined as expectation value for the upper limit of a given confidence interval under the assumption of no oscillations, is again obtained by simulations: $`\mathrm{sin}^2(2\mathrm{\Theta })=0.6(1.6,2.2)10^3`$ for 68 % C.L. (90% C.L., 95% C.L.), respectively. The error on the contour lines is given by the statistical error of the estimator’s cut value $`\mathrm{\Delta }\mathrm{ln}(\delta )`$ and therefore depends on the amount of simulated samples per hypothesis and the degree $`\delta `$ of confidence under consideration. For example, determining $`\mathrm{\Delta }\mathrm{ln}(90\%)`$ from a distribution of 1000 simulated samples relies on the upper 10 % tail or the highest 100 $`\mathrm{\Delta }\mathrm{ln}`$ values. We derive the error on $`\mathrm{\Delta }\mathrm{ln}`$ from the spread of the cut value taking the upper $`100\pm \sqrt{100}`$ samples. For KARMEN 2, the typical relative error on $`\mathrm{\Delta }\mathrm{ln}(90\%)`$ is in the range of $`23`$ %.

### III.2 LSND

For the LSND likelihood function, the method to deduce the confidence regions is identical to the one described in section III.1 for KARMEN 2. Technically, however, there are differences due to the large amount of computing time. This arises from the fact that the event samples are much larger (about 1000 events, depending on the oscillation parameters, instead of about 15) and the integration over different background contributions according to equation (6) is performed in as much as 255 discrete steps of different backgrounds BUB and BRB.

Therefore, we restricted the hypotheses tested from a fine grid for KARMEN 2 to some grid points representing the typical variation of expected oscillation events. For a fixed oscillation parameter $`\mathrm{\Delta }m^2=1`$eV2/c4, 4 oscillation scenarios were simulated with 1000 individual samples each (see Table 3).

As verified by earlier calculations NJP  and the KARMEN 2 analysis, the main variation of $`\mathrm{\Delta }\mathrm{ln}(\delta )`$ for a given confidence $`\delta `$ in the parameter space ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) follows the variation of the signal strength, i.e. the number of oscillation events $`N_{osc}`$ among the event sample. We therefore extrapolated the values of $`\mathrm{\Delta }\mathrm{ln}`$ in Table 3 linearly with $`N_{osc}`$ and adopted, at each grid point ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) with its expected number of oscillation events $`N_{osc}`$, the estimator value $`\mathrm{\Delta }\mathrm{ln}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)`$. Having constructed the estimator cut distribution this way, the last step is analogous to that used for KARMEN 2, i.e. to cut the LSND likelihood function with the given values $`\mathrm{\Delta }\mathrm{ln}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)`$ below the maximum and draw the contour lines. Figure 6 shows the contour lines of areas containing 68 %, 90 % and 95 % confidence.

Although parameter combinations with $`\mathrm{\Delta }m^2>10`$eV2/c4 are part of the 90% C.L. area, they are scarcely included in the more stringent interval of 68 % confidence, demonstrating the lower likelihood values for such oscillation scenarios.

It is useful to compare the extracted confidence regions with the favored regions given in LSNDfinal  where the likelihood function had been calculated based on a larger event sample (5697 events with $`20<E_e<200`$ MeV and no $`R_\gamma `$–cut). In addition, a simplified extraction of confidence regions according to equation (11) had been applied 666In fact, due in part to the larger statistical sample, in LSNDfinal  the results of such a simplified procedure were shown to be similar to the results obtained from an approximation to a full unified frequentist approach..

Figure 7 shows, that for lower values of $`\mathrm{\Delta }m^2`$, there is only a slight shift to larger oscillation amplitudes of the new analysis. A large overlap of the 90% C.L. area with the region according to $`\mathrm{\Delta }\mathrm{ln}=2.3`$ demonstrates the good statistical agreement of the analyses of both event samples. For higher values of $`\mathrm{\Delta }m^2`$, the differences are more pronounced. Compared to the older analysis, the confidence region is not only shifted considerably to larger amplitudes, but the likelihood values themselves are larger relative to the likelihood maximum which results in a larger area of 90% C.L.. These effects can be understood by the upper energy cut $`E_e<60`$ MeV applied in the actual analysis. Restricting the analysis mainly to the oscillation channel $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ favors also large values of $`\mathrm{\Delta }m^2`$ whereas a signal from $`\nu _\mu `$$``$$`\nu _e`$, derived in LSNDfinal  to be much smaller, but detectable mainly at higher energies, decreases the likelihood for such solutions as well as the absolute strength or oscillation amplitude.

Finally, we note the importance to derive values of $`\mathrm{\Delta }\mathrm{ln}`$ leading to correct coverage of confidence regions. These values differ considerably from the simplified Gaussian approach of a constant of 2.3 units. In an earlier analysis of the LSND data NJP , the typical estimator cut values were $`\mathrm{\Delta }\mathrm{ln}(90\%)2.53.5`$ on the basis of about 3000 events. Now, with even smaller event samples, the values of $`\mathrm{\Delta }\mathrm{ln}`$ further increase to $`3.54.5`$. This dependance of $`\mathrm{\Delta }\mathrm{ln}`$ is also underlined by the typical values for the much smaller KARMEN 2 sample size, where $`\mathrm{\Delta }\mathrm{ln}(90\%)`$ ranges from 3 to 5.

## IV Compatibility Analysis

Having analysed both experiments’ likelihood functions with the same consistent method, we can now use this method and its results to deduce quantitative statements on the question of statistical compatibility of both experimental results and, in the case of such compatibility, on the common parameter combinations ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$).

### IV.1 Level of compatibility

One of the most misleading but nevertheless very frequently used interpretation of the LSND and KARMEN results is to take the LSND region left over from the KARMEN exclusion curve as area of ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) favored by both experiments. Such an interpretation, though appealingly straight forward, completely ignores the information of both likelihood functions and reduces them to two discrete levels of a specific confidence $`\delta `$. To be able to correctly combine the two experimental results and extract the combined confidence regions, we have to use the original estimator distributions for KARMEN and for LSND. In a first step, the level of statistical compatibility of the two experimental outcomes has to be quantified. In the second step, the parameters favored by both experiments will be deduced.

In the following analysis, an experiment or event, in logical terms of frequency or probability of occurence, means repeating both, the KARMEN and the LSND experiment. Compatibility is then achieved if there are at least some oscillation parameters being the element of both the repeated KARMEN-like and LSND-like experiment’s confidence region.

We make the well justified assumption that the two experiments LSND and KARMEN are independent. Then, a two dimensional estimator distribution can be constructed for each hypothesis $`(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)_H`$ from the individual distributions $`\mathrm{\Delta }\mathrm{ln}_K`$ of KARMEN and $`\mathrm{\Delta }\mathrm{ln}_L`$ of LSND by an inverse projection. For each hypothesis, the pair of experimental values $`(\mathrm{\Delta }\mathrm{ln}_{KARMEN2},\mathrm{\Delta }\mathrm{ln}_{LSND})`$ is checked to be part of the ’inner’ $`\delta `$ % of the simulated 2–dimensional estimator distribution. There are different methods of ordering these 2–dimensional distributions all leading to very similar results NJP . In the following, we require each experimental value $`\mathrm{\Delta }\mathrm{ln}_{KARMEN2}`$ and $`\mathrm{\Delta }\mathrm{ln}_{LSND}`$ to be contained in the corresponding 1–dimensional interval of confidence. This selection criterion is equivalent to selecting a rectangle of frequency of the 2–dim estimator distribution. The resulting confidence is $`\delta =\delta _{KARMEN}\delta _{LSND}`$. Since we do not weight the experiments, we define $`\delta _{KARMEN}=\delta _{LSND}`$, so that for a combined confidence of e.g. 81 %, we demand $`\mathrm{\Delta }\mathrm{ln}_{KARMEN2}<\mathrm{\Delta }\mathrm{ln}_K(90\%)`$ and $`\mathrm{\Delta }\mathrm{ln}_{LSND}<\mathrm{\Delta }\mathrm{ln}_L(90\%)`$.

At a level of $`\delta _{KARMEN}=\delta _{LSND}=60\%`$, no parameter combination ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) fulfills the requirement $`(\mathrm{\Delta }\mathrm{ln}_{KARMEN2}<\mathrm{\Delta }\mathrm{ln}_K)(\mathrm{\Delta }\mathrm{ln}_{LSND}<\mathrm{\Delta }\mathrm{ln}_L)`$ any more. This corresponds to the fact that, at the level of combined confidence

$$\delta _{inc}=0.6^2=36\%,$$
(12)
the two experiments are completely incompatible. Coming back to our definition of probability of occurence, in 64 % of repetitions of double experiments, a KARMEN 2 plus an LSND experiment with their typical parameters and statistics, the outcome would consist of statistically compatible results with a specified set of oscillation parameters ($`\mathrm{sin}^2(2\mathrm{\Theta })`$, $`\mathrm{\Delta }m^2`$).

### IV.2 Common oscillation parameters

In the following, we analyse the preferred regions of oscillation parameters in the case of compatibility defined in the section before. If two experiments are independent producing two likelihood functions as result of the event analysis, the recipe for a combined analysis is straight forward. The likelihoods can be multiplied, i.e. in our case the logarithms of the likelihood functions are added. The absolute values of the likelihood functions are somewhat arbitrary. A recommended presentation of $`\mathrm{ln}`$ is to normalize the individual functions $`\mathrm{ln}_{KARMEN}`$ and $`\mathrm{ln}_{LSND}`$ to a point in ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) where they are equally sensitive to a potential signal dagos2 . In our case of the oscillation search this corresponds to values of $`\mathrm{sin}^2(2\mathrm{\Theta })=0`$. Therefore we define the combined likelihood function

$`\mathrm{ln}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)`$ $`=`$ $`\left\{\mathrm{ln}_{KARMEN2}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)\mathrm{ln}_{KARMEN2}(\mathrm{sin}^2(2\mathrm{\Theta })=0)\right\}`$ (13)
$`+`$ $`\left\{\mathrm{ln}_{LSND}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)\mathrm{ln}_{LSND}(\mathrm{sin}^2(2\mathrm{\Theta })=0)\right\}.`$
Figure 8 shows the combined function $`\mathrm{ln}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)`$ with a maximum of $`\mathrm{ln}`$($`\mathrm{sin}^2(2\mathrm{\Theta })`$$`=`$1,$`\mathrm{\Delta }m^2`$$`=`$0.05)$`=`$21.5 on a long flat ’ridge’ of low $`\mathrm{\Delta }m^2`$ values. The positive signal of LSND dominates the combined likelihood. However, the maximum in comparison to the no–oscillation value is reduced (see Equ. 8) as well as the shape at higher values of $`\mathrm{\Delta }m^2`$ is significantly modified due to the KARMEN 2 likelihood function.

To determine the confidence regions in ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$), again we apply the unified approach, i.e. we create the estimator distribution $`\mathrm{\Delta }\mathrm{ln}`$K+L for various hypotheses ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$)H and compare the experiment’s value $`\mathrm{\Delta }\mathrm{ln}`$KARMEN2+LSND with the cut value $`\mathrm{\Delta }\mathrm{ln}`$$`{}_{K+L}{}^{}(\delta )`$ for a given confidence. To get the estimator, a simulated event sample now consists of a KARMEN 2-like and an LSND-like MC sample each analysed with the appropriate likelihood definitions. Then both likelihood functions are added following equation (13) and the estimator $`\mathrm{\Delta }\mathrm{ln}_{K+L}=\mathrm{ln}_{K+L}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)_{max}\mathrm{ln}_{K+L}(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)_H`$ is calculated.

The typical cut values for a confidence of $`\delta =90`$ % are slightly higher than the individual LSND ones, ranging from 3.5 to 5.5. Figure 9 shows the confidence regions of the oscillation parameters for the combined likelihood analysis. The total confidence of a parameter region is hereby reduced by the fraction of incompatibility of the two experiments, so

$$\delta _{tot}=\delta (1\delta _{inc})$$
(14)
which results in e.g. $`\delta _{tot}=0.90.64=58`$ % total confidence for the parameter combinations within the area denoted as 90% C.L. area in Figure 9.

The combination of both experiments’ results with this consistent frequentist approach demonstrates that $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ solutions with $`\mathrm{\Delta }m^2>10`$eV2/c4 are excluded. There remain essentially two solutions in the parameter space of $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations, one at $`\mathrm{\Delta }m^27`$eV2/c4 and the area with $`\mathrm{\Delta }m^2<1`$eV2/c4. The latter one, though graphically rather large in Figure 9 corresponds to just one solution in the ’phase space’ of oscillations: Any $`\overline{\nu }_e`$ signal from oscillation has rather low energy, the evolution of the oscillation probability (equ. 1) being just at the beginning of the oscillation length. With an oscillation parameter $`\mathrm{\Delta }m^2=0.5`$eV2/c4, and a typical lower energy of $`E_\nu =20`$ MeV from $`\mu ^+`$ decay at rest, the first maximum of oscillations of $`\overline{\nu }_\mu `$ into $`\overline{\nu }_e`$ would be at a distance

$$d_{1.max}=\frac{\pi }{2}\frac{E_\nu }{1.27\mathrm{\Delta }m^2}50\mathrm{m}$$
(15)
from the neutrino source. For such oscillation parameters, the negative result of KARMEN 2 compared to the excess of LSND mainly reflects the different detector distances $`d_{KARMEN}17`$ m versus $`d_{LSND}30`$ m from the source. This is one way to reconcile the results of the two experiments.

To summarize, at a joint confidence level of $`\delta _{tot}=0.36+0.640.6880`$ %, the two experiments are either incompatible or lead to just one single oscillation solution at $`\mathrm{\Delta }m^2<1`$eV2/c4. The remaining 20 % confidence separate into an enlargement of the parameter region of low $`\mathrm{\Delta }m^2`$ values and into a second oscillation solution at $`\mathrm{\Delta }m^27`$eV2/c4.

## V Conclusion

We have analysed both the KARMEN 2 and the LSND final data with a maximum likelihood method using a similar event–based likelihood function. The likelihood functions as functions of the free parameters $`\mathrm{sin}^2(2\mathrm{\Theta })`$ and $`\mathrm{\Delta }m^2`$ demonstrate the different outcomes of the experiments, namely a clear, statistically significant excess of $`\overline{\nu }_e`$ in LSND versus no indication of an oscillation signal in KARMEN 2. To deduce regions of correct coverage, we applied a unified frequentist approach to both likelihood analyses individually. The results underline the feasibility of as well as the necessity for such an approach.

A quantitative joint statistical analysis has been performed leading to a level of 36 % incompatibility of the experimental outcomes, corresponding to individual confidence levels of 60 %. For the cases of statistical compatibility, the common parameter regions have been identified on the basis of the unified frequentist approach applied to the combined likelihood function of KARMEN 2 and LSND. The derived confidence regions in ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) clearly differ from an often applied but incorrect graphical overlap of the confidence regions of the individual experiments. There are two oscillation scenarios with either $`\mathrm{\Delta }m^27`$eV2/c4 or $`\mathrm{\Delta }m^2<1`$eV2/c4 compatible with both experiments.

We performed a joint statistical analysis incorporating some of the systematic uncertainties of the experiments, such as neutrino flux uncertainty, accuracy of known cross sections and resolution functions of both experiments.

Further –unknown– systematic uncertainties might become evident only by performing new experiments to confirm or discard the results described here. Figure 10 shows the intended sensitivity of a new experiment, MiniBooNE at Fermilab BooNE , which is under construction and will independently crosscheck the LSND evidence.

In addition, we did not incorporate the results of other experiments on the oscillation parameters investigated. As shown in Figure 10, the deduced favored region in ($`\mathrm{sin}^2(2\mathrm{\Theta })`$,$`\mathrm{\Delta }m^2`$) is partly overlaped by 90% C.L. exclusion curves of other experiments. In particular, the recent NOMAD $`\nu _\mu `$$``$$`\nu _e`$ result NOMAD  clearly reduces the overall likelihood for the $`\mathrm{\Delta }m^27`$eV2/c4 solution. A complete analysis should include these results on the basis of the same statistical method, i.e. a consistent frequentist analysis. This implies, however, the detailed knowledge of experimental data and resolution functions of these experiments not accessible to us. Furthermore, the exclusion curve from the Bugey experiment is based on the disappearance search $`\overline{\nu }_e`$$``$$`\overline{\nu }_x`$. Combining this experiment correctly with the appearance results of $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ or $`\nu _\mu `$$``$$`\nu _e`$ in terms of mixing angles would therefore also require a full three– or four–dimensional mixing scheme with theoretical models including sterile neutrinos bilenky ; Hirsch ; Ioann ; Ma ; Lorenz ; strumia  or CPT violation barenb  and would become model-dependent.

## VI Acknowledgements

We would like to thank the members of both collaborations for the encouragement and help to make this final joint analysis possible. In particular, we acknowledge the support of and many discussions with W.C. Louis concerning the LSND as well as G. Drexlin concerning the KARMEN experiment. Without the openness and trust to analyse both data sets in detail, such an analysis wouldn’t have been feasible. K.E. also wants to thank the Alexander von Humboldt foundation for financial support during visits to the Los Alamos National Laboratory. Clarifying discussions with Kai Zuber on the principles of combining experiments with contradictory central statements were very helpful.


format: text

text:
# Corrections to Scaling in the Droplet Picture of Spin Glasses

## Abstract

The energy of a droplet of linear extent $`l`$ in the “droplet theory” of spin glasses goes as $`l^\theta `$ for large $`l`$. It is argued that this formula needs to be modified by the addition of a scaling correction $`l^\omega `$ in order to accurately describe droplet energies at the length scales currently probed in numerical simulations. With this simple modification all equilibrium numerical data on Ising spin glasses in two, three and four dimensions becomes compatible with the droplet model.

The controversy as to the nature of the ordered state of spin glasses still remains to be resolved. There are two theories: the “droplet theory” McM ; BM ; FH  and the replica symmetry breaking (RSB) theory of Parisi Parisi ; sgb . Many numerical studies have been done in an attempt to resolve the controversy with results which are often uncomfortable for proponents of both theories. The most natural summary of the current numerical situation is the TNT (trivial, non-trivial) picture of Refs. KM ; PY . On this picture, which we shall explain in more detail below, the Parisi spin overlap function $`P(q,L)`$ for a system of linear dimension $`L`$ seems to behave as expected from the RSB theory i.e. in a non-trivial fashion, while the behavior of the link overlap seems to be more in accordance with droplet model ideas as its variance seems to be tending to zero as $`L`$ increases to infinity (which is regarded as trivial). In this paper we shall show that a simple and natural correction to scaling term added to the usual droplet energy expression can explain the origin of the TNT picture.

In the RSB theory there are low-energy excitations in which a finite fraction of the spins in the system are reversed but which only cost a finite amount of energy in the thermodynamic limit. The surface of these excitations is space-filling so that the fractal dimension of their surface $`d_s`$ is the same as the space dimension $`d`$ Marinari . In the droplet theory, the lowest energy excitations of linear spatial extent $`l`$ typically costs of order $`l^\theta `$, where $`\theta `$ is a positive exponent, of order 0.20 in three dimensions 3D  and 0.70 in four dimensions 4D  according to studies of the energy to create a domain wall right across the system. Thus in the thermodynamic limit the excitations which flip a finite fraction of the spins cost an infinite amount of energy and also the surface of these excitations is not space filling, as $`d_s<d`$.

The Hamiltonian which is usually studied in numerical work is the Ising spin glass model:

$$=\underset{<i,j>}{}J_{ij}S_iS_j,$$
(1)
where the sites $`i`$ lie on the sites of a cubic lattice with $`N=L^d`$ sites, $`S_i=\pm 1`$, the $`J_{ij}`$ have a Gaussian distribution of zero mean and unit variance and couple nearest-neighbor sites on the lattice. Typically for equilibrium studies at low temperatures (done at temperatures $`T`$ much less than the critical temperature $`T_c`$ in order to diminish the complications from critical point effects Drossel ), $`L14`$ for 3D and $`L5`$ for 4D. By studying two real replicas of the system, with spins $`S_i`$ and $`S_i^{^{}}`$ one can derive the Parisi spin overlap function $`P(q,L)`$, where

$$q=\frac{1}{N}\underset{i}{}S_iS_i^{^{}}$$
(2)
and study the properties of link overlaps

$$q_l=\frac{1}{N_b}\underset{<ij>}{}S_iS_jS_i^{^{}}S_j^{^{}},$$
(3)
where $`N_b`$ is the number of mearest-neighbor bonds on the lattice. On the droplet picture $`P(q,L)`$ should tend as $`L\mathrm{}`$ to a delta function at $`q=q_{EA}`$, and $`P(q,L)T/L^\theta ^{^{}}`$ for all $`qq_{EA}`$ where the expectation is that $`\theta ^{^{}}=\theta `$. In fact it is found that $`\theta ^{^{}}0`$ and that $`P(q,L)`$ looks very similar to what is predicted from the mean-field RSB picture Marinari  in that $`P(q,L)0`$ for all $`qq_{EA}`$ It is this result which advocates of the droplet model have found hardest to understand.

Studies of link overlaps are more in accordance with droplet model expectations. Their variance is predicted on the droplet model to decrease with $`L`$ as $`T/L^{\mu _l}`$ with $`\mu _l=2(dd_s)+\theta ^{^{}}`$ Drossel2 . The most recent investigations of link overlaps have been performed by adding a bulk perturbation $`ϵq_l`$ to $``$ and studying its effect in the ground state (see below). Using this method Palassini and Young PY  report that in 3D $`d_s=2.58\pm 0.02`$ and $`\theta ^{^{}}=0.02\pm 0.03`$ while in 4D $`d_s=3.77\pm 0.05`$ and $`\theta ^{^{}}=0.03\pm 0.05`$. The droplet model expectation is again that $`\theta ^{^{}}=\theta `$, which is clearly not satisfied at the $`L`$ values which can be currently studied. On the RSB picture the variance of $`q_l`$ should be non-zero for large L, which is not in accord with the results of PY . (However, Marinari and Parisi MP  did obtain a finite variance of the link overlap in the large $`L`$ limit by extrapolating it in an ad hoc way as $`1/L`$).

Thus if one accepts the TNT version of the numerical results but believes in the droplet model, the task is to understand why $`\theta ^{^{}}\theta `$ at least at the small values of $`L`$ which can be currently reached in equilibrium numerical studies. We will assume that the droplet picture is essentially correct and explain why corrections to scaling result in $`\theta ^{}`$ appearing to be close to zero for small $`L`$ values.

It is important to realise what the low-energy excitations look like in spin glasses. Consider a domain wall crossing a system $`M\times M\times L`$. The energy cost of this domain on the droplet picture will be of order $`M^\theta `$ and the domain wall will be fractal with dimension $`d_s`$ and with an area of order $`M^{d_s}`$. In three or more dimensions the interface may have holes throught it. Because it is fractal the extent of this wandering is of order $`M`$. The wandering of the interface by an amount of order $`M`$ affects the determination of $`\theta `$, see Ref. Carter , where it was found that the best results were obtained when $`LM`$ as then the interface is not affected by its interactions with the ends of the system. Suppose now we have two domain walls across the system. Then if their separation is large compared to $`M`$ they will be unaffected by the presence of the other. However, if they are closer together then they interfere with each other and their overall energy will be greater than if the other one were absent. In other words, domain walls effectively repel each other. If they have a separation of order $`l`$, then the repulsive energy between them would be expected to vary with their separation as a power law, $`l^\omega ^{^{}}`$. No investigations of $`\omega ^{^{}}`$ seem to exist in the literature.

Consider now a droplet of linear extent $`l`$. It too will have a fractal surface described by $`d_s`$. Pictures of large droplets have now appeared KM , and a systematic investigation of them is in HKM . Their fractal nature ensures that in three or more dimensions that they have holes through them, giving them a sponge-like appearance. However, because the surface of the droplet may wander by a distance of order $`l`$ the energy of the droplet $`E`$ will be modified by its wandering and “collisions” with itself to a form which we suppose by analogy with the above is

$$E=Al^\theta +Bl^\omega .$$
(4)
The term $`Bl^\omega `$ is a scaling correction to the form of the droplet energy at large $`l`$ and $`A`$ and $`B`$ are positive constants. We shall argue that with this scaling correction it is possible to understand the TNT description of the numerical scene. The exponent $`\omega `$ is only well-defined as a correction to scaling exponent when $`l\mathrm{}`$, but we shall assume that Eq. (4) has utility outside this limit. We do not know whether $`\omega ^{^{}}=\omega `$, but fortunately there already exists numerical data on the value of $`\omega `$ in 3D from the work of Lamarcq et al. LQ .

These authors first computed the ground state for both $`N=6^3`$ and $`N=10^3`$. They then chose an arbitrary reference spin and flipped it along with a cluster containing $`v1`$ other spins connected to it. They next minimized the energy of this cluster by exchange Monte Carlo, but with the constraint that the reference spin is held fixed and the cluster was always connected and of size $`v`$. They found the largest extension (mean end-to-end distance) $`l`$ of the cluster and found that for $`v33`$ its energy $`E(l)`$ varied as $`l^\omega `$, where the exponent $`\omega `$ was $`0.13\pm 0.02`$. For larger $`v`$ values they found that $`E(l)`$ was increasing rather than decreasing with $`l`$ but they were unsure whether this might not be an artifact of insufficient numbers of Monte Carlo steps. Note that according to Eq. (4), when both $`\theta `$ and $`\omega `$ are small $`E(l)`$ will have a shallow minimum. The clusters generated by their procedure were not compact (that is, $`vl^d`$), but on the droplet picture, which is a scaling picture associated with the zero-temperature fixed point McM ; BM  it is only required that large clusters be compact.

Eq. (4) is more illuminating when expressed in variables associated with its minimum at $`l=l_0`$:

$$E(l)=C\left(\frac{(l/l_0)^\theta }{\theta }+\frac{(l_0/l)^\omega }{\omega }\right),$$
(5)
where $`C`$ is a positive constant. The position of the minimum is essentially unknown but there is a hint from the upturn in $`E(l)`$ seen in Ref. LQ  that $`l_0`$ might lie between 7 and 10 lattice spacings in 3D. Then because the values of both $`\theta `$ and $`\omega `$ are small there is only a weak dependence of $`E(l)`$ on $`l`$ in the region accessible to numerical studies (say 14 lattice spacings). The minimum is very shallow so that the apparent value of $`\theta `$, which is $`\theta ^{^{}}`$, will be close to zero. For example, the ratio $`E(14)/E(7)=1.006`$ if $`l_0`$ is 7, whereas in the absence of the correction to scaling term this ratio is 1.15. In 4D the values of $`\omega `$ and $`l_0`$ are undetermined at present, but provided $`l_0`$ happens to lie in the region open to numerically studies, then again $`\theta ^{^{}}`$ would be small.

While it is inevitable that exponents in 4D will be difficult to determine in numerical studies this is not so for 2D. There seems to be indirect evidence which supports an $`E(l)`$ as in Eq. (4) also for 2D systems. In two dimensions $`\theta `$ is very accurately determined as systems of $`480^2`$ can be studied HY . The exponent associated with a single domain wall, that is $`\theta `$, is $`0.282(2)`$ Carter ; HY . However, studies where, say, $`\theta `$ is determined from the effects of thermally excited droplets (such as Monte Carlo simulations of the spin-glass susceptibility) yield an apparent $`\theta ^{^{}}`$ close to $`0.47`$ KA . This discrepancy has long been a puzzle, and has prompted suggestions that perhaps different exponents describe domain wall energies and droplet energies KA . However, Eq. (4) with $`\omega 0.47`$ would seem to offer another way of resolving the descrepancy. For small values of $`l`$, droplets would apparently have energies decreasing as $`l^\omega `$ as their energies are dominated by the correction to scaling term, but at large values of $`l`$ the decrease with $`l`$ will be slower and be as $`l^{|\theta |}`$, as expected from the conventional droplet approach.

This crossover between the large and small $`l`$ dependencies can be seen in the work of Middleton M . He studied link overlaps as in PY  and MP , by comparing the unperturbed $`J_{ij}`$ ground state, $`\{S_i^0\}`$, with the $`ϵ`$-perturbed state, where $`J_{ij}J_{ij}ϵS_i^0S_j^0/N_b`$. The fraction of link values $`S_iS_j`$ on nearest-neighbor bonds $`<ij>`$ which are changed by the perturbation is $`(1q_l)`$. On the droplet model the sample average of this, $`\overline{1q_l}`$, is predicted to decrease as $`ϵ/L^{\mu _l}`$. Using the values $`\theta =0.28`$ and $`d_s=1.27`$, $`\mu _l1.18`$, but with $`\theta ^{^{}}`$ at its small $`l`$ value of $`0.47`$ the expected value of $`\mu _l0.99`$. Middleton studied the effective value of $`\mu _l`$ at length scale L using the definition

$$\mu _l^{\mathrm{eff}}=\frac{d[\mathrm{ln}(\overline{1q_l})]}{d\mathrm{ln}(L)}.$$
(6)
At $`L=16`$ $`\mu _l^{eff}`$ was found to be 0.99 and increased to the large $`l`$ expected value 1.18 when $`L200`$. In other words, $`\mu _l^{eff}`$ seems to interpolate beteen the two limits as expected on the basis of our correction to scaling.

Middleton himself attributed this gradual evolution of $`\mu _l^{eff}`$ with $`L`$ to finite size corrections which can arise when studying quantities which have contributions from droplets of all sizes up to the system size $`L`$. The domain wall energy exponent $`\theta `$ is free of this problem as by definition domain walls only exist on the scale $`L`$. He showed that even if $`E(l)`$ varies just as $`l^\theta `$ for all $`l`$, then the contributions of droplets on all length scales up to $`L`$ means that $`\mu _l^{eff}`$ approaches its asymptotic value only as $`1/L^{d\mu _l}`$. These scale averaging corrections decrease in 2D with the exponent $`d\mu _l0.82`$ and are therefore unlikely to be important when $`L100`$. In three dimensions $`d\mu _l1.3`$. As in 3D only modest values of $`L`$ can be studied the scale averaging corrections may there be a significant effect.

The effect of the correction to scaling in the formula for the droplet energy will have relevance beyond the confines of numerical simulations. The droplet model provides scaling laws for time-dependent quantities such as dynamical susceptibilities in terms of a length scale $`L(t)`$ which increases in a logarithmic manner with time $`t`$ due to thermal activation of droplets on the length scale $`L(t)`$. In practice the length scale which can be explored in real experiments is rather limited (typically more than 10 lattice spacings but probably less than 200 JYN ). The fact that this length scale is not very large means that corrections to asymptotic formula will be important and that will include the correction to the scaling energy proposed in this paper.

To summarize: there exists a simple correction to the usual scaling formula for the energy of a droplet which is consistent with the TNT summary of the numerical data on spin glasses. It would be useful to have direct studies of the two correction to scaling exponents exponents $`\omega `$ and $`\omega ^{^{}}`$. In 3D and 4D, information on the value of $`l`$ at the minimum of $`E(l)`$, $`l_0`$, is needed. As perhaps $`l_0`$ might have been already seen in the 3D work of Ref. LQ , a precise determination could be possible with a little extra effort. In 2D an accurate determination of both exponents might be achievable with the techniques used in Refs. HY  and M .

###### Acknowledgements.

I would like to thank Alan Bray for many discussions.


text:
# Mathematical Structure of Rabi Oscillations in the Strong Coupling Regime

## 1 Introduction

Coherent states or generalized coherent states play an important role in quantum physics, in particular, quantum optics, see  and . They also play an important one in mathematical physics. See the textbook . For example, they are very useful in performing stationary phase approximations to path integral, , , .

Coherent operators which produce coherent states are very useful because they are unitary and easy to handle. The basic reason is probably that they are subject to the elementary Baker-Campbell-Hausdorff (BCH) formula. Many basic properties of them are well–known, see  or .

Generalized coherent operators which produce generalized coherent states are also useful. But they are not so easy to handle in spite of having the disentangling one corresponding to the elementary BCH formula. In  and  the author determined all matrix elements of generalized coherent operators based on Lie algebras su(1,1) and su(2). They are interesting by themselves, but moreover have a very interesting application.

In  Frasca dealt with the Jaynes–Cummings model which describes a two–level atom interacting with a single radiation mode (see  for a general review) in the strong coupling regime (not weak coupling one !) and showed that Rabi frequencies are obtained by matrix elements of coherent operator under the rotating–wave approximation. His aim was to explain the recent experimental finding on Josephson junctions .

This is an interesting result and moreover his method can be widely generalized. See also  for an another example dealt with in the strong coupling regime.

In this paper we generalize the Jaynes–Cummings Hamiltonian by making use of some operators based on Lie algebras su(1,1) and su(2), and study a mathematical structure of Rabi floppings of these extended models in the strong coupling regime.

We show that (generalized) Rabi frequencies are also given by matrix elements of generalized coherent operators under the rotating–wave approximation. We believe that the results will give a new aspect to Quantum Optics or Mathematical Physics.

We also apply our results to the construction of controlled unitary gates in Quantum Computation in the last section.

Lastly we discuss an application to Holonomic Quantum Computation, but our discussion is not complete.

## 2 Coherent and Generalized Coherent Operators

### 2.1 Coherent Operator

Let $`a(a^{})`$ be the annihilation (creation) operator of the harmonic oscillator. If we set $`Na^{}a`$ (: number operator), then

$$[N,a^{}]=a^{},[N,a]=a,[a^{},a]=\mathrm{𝟏}.$$
(1)
Let $``$ be a Fock space generated by $`a`$ and $`a^{}`$, and $`\{|n|n𝐍\{0\}\}`$ be its basis. The actions of $`a`$ and $`a^{}`$ on $``$ are given by

$$a|n=\sqrt{n}|n1,a^{}|n=\sqrt{n+1}|n+1,N|n=n|n$$
(2)
where $`|0`$ is a normalized vacuum ($`a|0=0\mathrm{and}0|0=1`$). From (2) state $`|n`$ for $`n1`$ are given by

$$|n=\frac{(a^{})^n}{\sqrt{n!}}|0.$$
(3)
These states satisfy the orthogonality and completeness conditions

$$m|n=\delta _{mn},\underset{n=0}{\overset{\mathrm{}}{}}|nn|=\mathrm{𝟏}.$$
(4)
Definition We call a state

$$|z=\text{e}^{za^{}\overline{z}a}|0U(z)|0\text{for}z𝐂$$
(5)
the coherent state.

### 2.2 Generalized Coherent Operator Based on $`su(1,1)`$

Let us state generalized coherent operators and states based on $`su(1,1)`$.

We consider a spin $`K(>0)`$ representation of $`su(1,1)sl(2,𝐂)`$ and set its generators $`\{K_+,K_{},K_3\}((K_+)^{}=K_{})`$,

$$[K_3,K_+]=K_+,[K_3,K_{}]=K_{},[K_+,K_{}]=2K_3.$$
(6)
We note that this (unitary) representation is necessarily infinite dimensional. The Fock space on which $`\{K_+,K_{},K_3\}`$ act is $`_K\{|K,n|n𝐍\{0\}\}`$ and whose actions are

$`K_+|K,n`$ $`=`$ $`\sqrt{(n+1)(2K+n)}|K,n+1,K_{}|K,n=\sqrt{n(2K+n1)}|K,n1,`$
$`K_3|K,n`$ $`=`$ $`(K+n)|K,n,`$ (7)
where $`|K,0`$ is a normalized vacuum ($`K_{}|K,0=0`$ and $`K,0|K,0=1`$). We have written $`|K,0`$ instead of $`|0`$ to emphasize the spin $`K`$ representation, see . From (2.2), states $`|K,n`$ are given by

$$|K,n=\frac{(K_+)^n}{\sqrt{n!(2K)_n}}|K,0,$$
(8)
where $`(a)_n`$ is the Pochammer’s notation $`(a)_na(a+1)\mathrm{}(a+n1).`$ These states satisfy the orthogonality and completeness conditions

$$K,m|K,n=\delta _{mn},\underset{n=0}{\overset{\mathrm{}}{}}|K,nK,n|=\mathrm{𝟏}_K.$$
(9)
Now let us consider a generalized version of coherent states :

Definition We call a state

$$|z=V(z)|K,0\text{e}^{zK_+\overline{z}K_{}}|K,0\text{for}z𝐂.$$
(10)
the generalized coherent state (or the coherent state of Perelomov’s type based on $`su(1,1)`$ in our terminology).

Here let us construct an example of this representation. First we set

$$K_+\frac{1}{2}\left(a^{}\right)^2,K_{}\frac{1}{2}a^2,K_3\frac{1}{2}\left(a^{}a+\frac{1}{2}\right),$$
(11)
then it is easy to check that these satisfy the commutation relations (6). That is, the set $`\{K_+,K_{},K_3\}`$ gives a unitary representation of $`su(1,1)`$with spin $`K=1/4\text{and}3/4`$, . Now we also call an operator

$$S(z)=\text{e}^{\frac{1}{2}\{z(a^{})^2\overline{z}a^2\}}\text{for}z𝐂$$
(12)
the squeezed operator, see the book .

### 2.3 Generalized Coherent Operator Based on $`su(2)`$

Let us state generalized coherent operators and states based on $`su(2)`$.

We consider a spin $`J(>0)`$ representation of $`su(2)sl(2,𝐂)`$ and set its generators $`\{J_+,J_{},J_3\}((J_+)^{}=J_{})`$,

$$[J_3,J_+]=J_+,[J_3,J_{}]=J_{},[J_+,J_{}]=2J_3.$$
(13)
We note that this (unitary) representation is necessarily finite dimensional. The Fock space on which $`\{J_+,J_{},J_3\}`$ act is $`_J\{|J,n|0n2J\}`$ and whose actions are

$`J_+|J,n`$ $`=`$ $`\sqrt{(n+1)(2Jn)}|J,n+1,J_{}|J,n=\sqrt{n(2Jn+1)}|J,n1,`$
$`J_3|J,n`$ $`=`$ $`(J+n)|J,n,`$ (14)
where $`|J,0`$ is a normalized vacuum ($`J_{}|J,0=0`$ and $`J,0|J,0=1`$). We have written $`|J,0`$ instead of $`|0`$ to emphasize the spin $`J`$ representation, see . From (2.3), states $`|J,n`$ are given by

$$|J,n=\frac{(J_+)^n}{\sqrt{n!{}_{2J}{}^{}P_{n}^{}}}|J,0.$$
(15)
These states satisfy the orthogonality and completeness conditions

$$J,m|J,n=\delta _{mn},\underset{n=0}{\overset{2J}{}}|J,nJ,n|=\mathrm{𝟏}_J.$$
(16)
Now let us consider a generalized version of coherent states :

Definition We call a state

$$|z=W(z)|J,0\text{e}^{zJ_+\overline{z}J_{}}|J,0\text{for}z𝐂.$$
(17)
the generalized coherent state (or the coherent state of Perelomov’s type based on $`su(2)`$ in our terminology).

A comment is in order. We can construct the spin $`K`$ and $`J`$ representations by making use of Schwinger’s boson method. But we don’t repeat here, see for example .

## 3 Matrix Elements of Coherent and Generalized Coherent Operators $`\mathrm{}`$

### 3.1 Matrix Elements of Coherent Operator

We list matrix elements of coherent operators $`U(z)`$.

The Matrix Elements The matrix elements of $`U(z)`$ are :

$`(\text{i})nmn|U(z)|m=\text{e}^{\frac{1}{2}|z|^2}\sqrt{{\displaystyle \frac{n!}{m!}}}(\overline{z})^{mn}L_{n}^{}{}_{}{}^{(mn)}(|z|^2),`$ (18)
$`(\text{ii})nmn|U(z)|m=\text{e}^{\frac{1}{2}|z|^2}\sqrt{{\displaystyle \frac{m!}{n!}}}z^{nm}L_{m}^{}{}_{}{}^{(nm)}(|z|^2),`$ (19)
where $`L_{n}^{}{}_{}{}^{(\alpha )}`$ is the associated Laguerre’s polynomial defined by

$$L_{k}^{}{}_{}{}^{(\alpha )}(x)=\underset{j=0}{\overset{k}{}}(1)^j\left(\genfrac{}{}{0pt}{}{k+\alpha }{kj}\right)\frac{x^j}{j!}.$$
(20)
In particular $`L_kL_{k}^{}{}_{}{}^{(0)}`$ is the usual Laguerre’s polynomial and these are related to diagonal elements of $`U(z)`$.

### 3.2 Matrix Elements of Coherent Operator Based on $`su(1,1)`$

We list matrix elements of $`V(z)`$ coherent operators based on $`su(1,1)`$. In this case it is always $`2K>1`$ ($`2K=1`$ under some regularization).

The Matrix Elements The matrix elements of $`V(z)`$ are :

$`(\text{i})nm`$ $`K,n|V(z)|K,m=\sqrt{{\displaystyle \frac{n!m!}{(2K)_n(2K)_m}}}(\overline{\kappa })^{mn}(1+|\kappa |^2)^{K\frac{n+m}{2}}\times `$ (21)
$`{\displaystyle \underset{j=0}{\overset{n}{}}}(1)^{nj}{\displaystyle \frac{\mathrm{\Gamma }(2K+m+nj)}{\mathrm{\Gamma }(2K)(mj)!(nj)!j!}}(1+|\kappa |^2)^j(|\kappa |^2)^{nj},`$
$`(\text{ii})nm`$ $`K,n|V(z)|K,m=\sqrt{{\displaystyle \frac{n!m!}{(2K)_n(2K)_m}}}\kappa ^{nm}(1+|\kappa |^2)^{K\frac{n+m}{2}}\times `$ (22)
$`{\displaystyle \underset{j=0}{\overset{m}{}}}(1)^{mj}{\displaystyle \frac{\mathrm{\Gamma }(2K+m+nj)}{\mathrm{\Gamma }(2K)(mj)!(nj)!j!}}(1+|\kappa |^2)^j(|\kappa |^2)^{mj},`$
where

$$\kappa \frac{\text{sinh}(|z|)}{|z|}z=\text{cosh}(|z|)\zeta .$$
(23)
The author doesn’t know whether or not the right hand sides of (21) and (22) could be written by making use of some special functions such as generalized Laguerre’s functions in (20). Therefore we set temporarily

$$F_{m}^{}{}_{}{}^{(nm)}(x:2K)=\underset{j=0}{\overset{m}{}}(1)^{mj}\frac{\mathrm{\Gamma }(2K+m+nj)}{\mathrm{\Gamma }(2K)(mj)!(nj)!j!}(1+x)^jx^{mj}$$
(24)
and $`F_{m}^{}{}_{}{}^{(0)}(x;2K)=F_m(x;2K)`$.

### 3.3 Matrix Elements of Coherent Operator Based on $`su(2)`$

We list matrix elements of $`W(z)`$ coherent operators based on $`su(2)`$. In this case it is always $`2J𝐍`$.

Matrix Elements The matrix elements of $`W(z)`$ are :

$`(\text{i})nm`$ $`J,n|W(z)|J,m=\sqrt{{\displaystyle \frac{n!m!}{{}_{2J}{}^{}P_{n}^{}{}_{2J}{}^{}P_{m}^{}}}}(\overline{\kappa })^{mn}(1|\kappa |^2)^{J\frac{n+m}{2}}\times `$ (25)
$`{\displaystyle \underset{j=0}{\overset{n}{}}}{}_{}{}^{}(1)_{}^{nj}{\displaystyle \frac{(2J)!}{(2Jmn+j)!(mj)!(nj)!j!}}(1|\kappa |^2)^j(|\kappa |^2)^{nj},`$
$`(\text{ii})nm`$ $`J,n|W(z)|J,m=\sqrt{{\displaystyle \frac{n!m!}{{}_{2J}{}^{}P_{n}^{}{}_{2J}{}^{}P_{m}^{}}}}\kappa ^{nm}(1|\kappa |^2)^{J\frac{n+m}{2}}\times `$ (26)
$`{\displaystyle \underset{j=0}{\overset{m}{}}}{}_{}{}^{}(1)_{}^{mj}{\displaystyle \frac{(2J)!}{(2Jmn+j)!(mj)!(nj)!j!}}(1|\kappa |^2)^j(|\kappa |^2)^{mj},`$
where

$$\kappa \frac{\text{sin}(|z|)}{|z|}z=\text{cos}(|z|)\eta .$$
(27)
Here $`_{}`$ means a summation over $`j`$ satisfying $`2Jmn+j0`$.

The author doesn’t know whether or not the right hand sides of (25) and (26) could be written by making use of some special functions. We set temporarily

$$F_{m}^{}{}_{}{}^{(nm)}(x:2J)=\underset{j=0}{\overset{m}{}}{}_{}{}^{}(1)_{}^{mj}\frac{(2J)!}{(2Jmn+j)!(mj)!(nj)!j!}(1x)^jx^{mj}$$
(28)
and $`F_{m}^{}{}_{}{}^{(0)}(x;2J)=F_m(x;2J)`$.

## 4 Jaynes–Cummings Models in the Strong Coupling Regime

In  Frasca treated the Jaynes–Cummings model and developped some method to calculate Rabi frequencies in the strong coupling regime. We in this section generalize the model and method, and show that Rabi frequencies in our extended model are given by matrix elements of generalized coherent operators under the rotating–wave approximation. This gives a unified approach to them.

Let $`\{\sigma _1,\sigma _2,\sigma _3\}`$ be Pauli matrices and $`\mathrm{𝟏}_2`$ a unit matrix :

$$\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\mathrm{𝟏}_2=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(29)
The Hamiltonian adopted in  is

$$\text{(N)}H_N=\omega \mathrm{𝟏}_2a^{}a+\frac{\mathrm{\Delta }}{2}\sigma _3\mathrm{𝟏}+g\sigma _1(a^{}+a)$$
(30)
where $`\omega `$ is the frequency of the radiation mode, $`\mathrm{\Delta }`$ the separation between the two levels of the atom, $`g`$ the coupling between the radiation field and the atom.

Moreover we want to treat the following Hamiltonians (our extension)

$`\text{(K)}H_K`$ $`=`$ $`\omega \mathrm{𝟏}_2K_3+{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma _3\mathrm{𝟏}_K+g\sigma _1(K_++K_{}),`$ (31)
$`\text{(J)}H_J`$ $`=`$ $`\omega \mathrm{𝟏}_2J_3+{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma _3\mathrm{𝟏}_J+g\sigma _1(J_++J_{}).`$ (32)
To treat these three cases at the same time we set

$$\{L_+,L_{},L_3\}=\{\begin{array}{cc}\text{(N)}\{a^{},a,N\}\hfill & \\ \text{(K)}\{K_+,K_{},K_3\}\hfill & \\ \text{(J)}\{J_+,J_{},J_3\}\hfill & \end{array}$$
(33)
and

$$H=H_0+V=\omega \mathrm{𝟏}_2L_3+\frac{\mathrm{\Delta }}{2}\sigma _3\mathrm{𝟏}_L+g\sigma _1(L_++L_{})$$
(34)
where we have written $`H`$ instead of $`H_L`$ for simplicity.

Mysteriously enough we cannot solve these simple models completely (maybe non–integrable), nevertheless we have found these models have a very rich structure.

For these (non–integrable) models we usually have two perturbation approaches :

Weak Coupling Regime ($`0<g\mathrm{\Delta }`$)

$$H_0=\omega \mathrm{𝟏}_2L_3+\frac{\mathrm{\Delta }}{2}\sigma _3\mathrm{𝟏}_L,V=g\sigma _1(L_++L_{}).$$
(35)
Strong Coupling Regime ($`0<\mathrm{\Delta }g`$)

$$H_0=\omega \mathrm{𝟏}_2L_3+g\sigma _1(L_++L_{}),V=\frac{\mathrm{\Delta }}{2}\sigma _3\mathrm{𝟏}_L.$$
(36)
In the following we consider only the strong coupling regime (see  for the weak one). First let us solve $`H_0`$ which is a relatively easy task.

Let $`W`$ be a Walsh–Hadamard matrix

$$W=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=W^1$$

then we can diagonalize $`\sigma _1`$ by using this $`H`$ as $`\sigma _1=W\sigma _3W^1`$. The eigenvalues of $`\sigma _1`$ is $`\{1,1\}`$ with eigenvectors

$$|1=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ 1\end{array}\right),|1=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ 1\end{array}\right)|\lambda =\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ \lambda \end{array}\right).$$

We note that

$`|11|={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=W\left(\begin{array}{cc}1& \\ & 0\end{array}\right)W^1,`$ (46)
$`|11|={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)=W\left(\begin{array}{cc}0& \\ & 1\end{array}\right)W^1,`$
$``$ $`|\lambda \lambda |={\displaystyle \frac{1}{2}}\left(\begin{array}{cc}1& \lambda \\ \lambda & 1\end{array}\right)=W\left(\begin{array}{cc}\frac{1+\lambda }{2}& \\ & \frac{1\lambda }{2}\end{array}\right)W^1.`$ (51)
Then we have

$`H_0`$ $`=`$ $`(W\mathrm{𝟏}_L)\left(\omega \mathrm{𝟏}_2L_3+g\sigma _3(L_++L_{})\right)(W^1\mathrm{𝟏}_L)`$
$`=`$ $`(W\mathrm{𝟏}_L)\left(\begin{array}{cc}\omega L_3+g(L_++L_{})& \\ & \omega L_3g(L_++L_{})\end{array}\right)(W^1\mathrm{𝟏}_L).`$
$`=`$ $`|11|\{\omega L_3+g(L_++L_{})\}+|11|\{\omega L_3g(L_++L_{})\}`$
$`=`$ $`{\displaystyle \underset{\lambda }{}}|\lambda \lambda |\{\omega L_3+\lambda g(L_++L_{})\}`$
$`=`$ $`{\displaystyle \underset{\lambda }{}}|\lambda \lambda |\left\{\text{e}^{\frac{\lambda x}{2}(L_+L_{})}\left(\mathrm{\Omega }L_3\right)\text{e}^{\frac{\lambda x}{2}(L_+L_{})}\right\}`$
$`=`$ $`{\displaystyle \underset{\lambda }{}}\left(|\lambda \text{e}^{\frac{\lambda x}{2}(L_+L_{})}\right)\left(\mathrm{\Omega }L_3\right)\left(\lambda |\text{e}^{\frac{\lambda x}{2}(L_+L_{})}\right)`$ (55)
where we have used the following

Key Formulas For $`\lambda =\pm 1`$ we have

$`(N)\omega a^{}a+\lambda g(a^{}+a)=\mathrm{\Omega }\text{e}^{\frac{\lambda x}{2}(a^{}a)}\left(N{\displaystyle \frac{g^2}{\omega ^2}}\right)\text{e}^{\frac{\lambda x}{2}(a^{}a)}`$
$`\text{where}\mathrm{\Omega }=\omega ,x=2g/\omega ,`$ (56)
$`(K)\omega K_3+\lambda g(K_++K_{})=\mathrm{\Omega }\text{e}^{\frac{\lambda x}{2}(K_+K_{})}K_3\text{e}^{\frac{\lambda x}{2}(K_+K_{})}`$
$`\text{where}\mathrm{\Omega }=\omega \sqrt{1(2g/\omega )^2},x=\text{tanh}^1(2g/\omega ),`$ (57)
$`(J)\omega J_3+\lambda g(J_++J_{})=\mathrm{\Omega }\text{e}^{\frac{\lambda x}{2}(J_+J_{})}J_3\text{e}^{\frac{\lambda x}{2}(J_+J_{})}`$
$`\text{where}\mathrm{\Omega }=\omega \sqrt{1+(2g/\omega )^2},x=\text{tan}^1(2g/\omega ).`$ (58)
The proof is not difficult, so we leave it to the readers. That is, we could diagonalize the Hamiltonian $`H_0`$. This is two–fold degenerate and its eigenvalues and eigenvectors are given respectively

$$(\text{Eigenvalues},\text{Eigenvectors})=\{\begin{array}{cc}(N)\omega n\frac{g^2}{\omega },|\lambda \text{e}^{\frac{\lambda x}{2}(a^{}a)}|n\hfill & \\ (K)\mathrm{\Omega }(K+n),|\lambda \text{e}^{\frac{\lambda x}{2}(K_+K_{})}|K,n\hfill & \\ (J)\mathrm{\Omega }(J+n),|\lambda \text{e}^{\frac{\lambda x}{2}(J_+J_{})}|J,n\hfill & \end{array}$$
(59)
for $`\lambda =\pm 1`$ and $`n𝐍\{0\}`$. For the latter convenience we set

$$\text{Eigenvalues}=\{E_n\},\text{Eigenvectors}=\{|\{\lambda ,n\}\}.$$
(60)
Then (4) can be written as

$$H_0=\underset{\lambda }{}\underset{n}{}E_n|\{\lambda ,n\}\{\lambda ,n\}|.$$
(61)
Next we would like to solve the following Schrödinger equation :

$$i\frac{d}{dt}\mathrm{\Psi }=H\mathrm{\Psi }=\left(H_0+\frac{\mathrm{\Delta }}{2}\sigma _3\mathrm{𝟏}_L\right)\mathrm{\Psi },$$
(62)
where we have set $`\mathrm{}=1`$ for simplicity. To solve this equation we appeal to the method of constant variation. First let us solve

$$i\frac{d}{dt}\mathrm{\Psi }=H_0\mathrm{\Psi },$$
(63)
which general solution is given by

$$\mathrm{\Psi }(t)=U_0(t)\mathrm{\Psi }_0=\text{e}^{itH_0}\mathrm{\Psi }_0$$
(64)
where $`\mathrm{\Psi }_0`$ is a constant state. It is easy to see from (61)

$$U_0(t)=\text{e}^{itH_0}=\underset{\lambda }{}\underset{n}{}\text{e}^{itE_n}|\{\lambda ,n\}\{\lambda ,n\}|.$$
(65)
The method of constant variation goes as follows. Changing like $`\mathrm{\Psi }_0\mathrm{\Psi }_0(t),`$ we insert (64) into (230). After some algebra we obtain

$$i\frac{d}{dt}\mathrm{\Psi }_0=\frac{\mathrm{\Delta }}{2}U_{0}^{}{}_{}{}^{}(\sigma _3\mathrm{𝟏}_L)U_0\mathrm{\Psi }_0.$$
(66)
We have only to solve this equation. If we set

$$H_F=\frac{\mathrm{\Delta }}{2}U_{0}^{}{}_{}{}^{}(\sigma _3\mathrm{𝟏}_L)U_0,$$
(67)
then we have easily from (65)

$`H_F`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\lambda ,\mu }{}}{\displaystyle \underset{m,n}{}}\text{e}^{it(E_mE_n)}\{\lambda ,m\}|(\sigma _3\mathrm{𝟏}_L)|\{\mu ,n\}|\{\lambda ,m\}\{\mu ,n\}|`$ (68)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{m,n}{}}\text{e}^{it\mathrm{\Omega }(mn)}m|\text{e}^{\lambda x(L_+L_{})}|n|\{\lambda ,m\}\{\lambda ,n\}|`$
where we have used the relation $`\lambda |\sigma _3=\lambda |`$. Remind that $`|n`$ is respectively

$$|n=\{\begin{array}{cc}(N)|n\hfill & \\ (K)|K,n\hfill & \\ (J)|J,n\hfill & \end{array}$$

In this stage we meet matrix elements of the coherent and generalized coherent operators $`\text{e}^{\lambda x(L_+L_{})}`$ in section 3 ($`z=\overline{z}=\lambda x`$).

Here we divide $`H_F`$ into two parts

$$H_F=H_{F}^{}{}_{}{}^{^{}}+H_{F}^{}{}_{}{}^{^{\prime \prime }}$$

where

$`H_{F}^{}{}_{}{}^{^{}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{n}{}}n|\text{e}^{\lambda x(L_+L_{})}|n|\{\lambda ,n\}\{\lambda ,n\}|,`$ (69)
$`H_{F}^{}{}_{}{}^{^{\prime \prime }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{\stackrel{m,n}{mn}}{}}\text{e}^{it\mathrm{\Omega }(mn)}m|\text{e}^{\lambda x(L_+L_{})}|n|\{\lambda ,m\}\{\lambda ,n\}|.`$ (70)
Noting

$$n|\text{e}^{x(L_+L_{})}|n=n|\text{e}^{x(L_+L_{})}|n$$

by the results in section 3, $`H_{F}^{}{}_{}{}^{^{}}`$ can be written as

$$H_{F}^{}{}_{}{}^{^{}}=\frac{\mathrm{\Delta }}{2}\underset{n}{}n|\text{e}^{x(L_+L_{})}|n\left\{|\{1,n\}\{1,n\}|+|\{1,n\}\{1,n\}|\right\},$$

from which we can diagonalize $`H_{F}^{}{}_{}{}^{^{}}`$ as

$$H_{F}^{}{}_{}{}^{^{}}=\frac{\mathrm{\Delta }}{2}\underset{n}{}\underset{\sigma }{}n|\text{e}^{x(L_+L_{})}|n\sigma |\{\sigma ,\psi _n\}\{\sigma ,\psi _n\}|$$
(71)
if we define a new basis

$$|\{\sigma ,\psi _n\}=\frac{1}{\sqrt{2}}(\sigma |\{1,n\}+|\{1,n\}),\sigma =\pm 1.$$

These states can be seen as so–called Schrödinger cat states, . From these we have

$`|\{1,n\}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\{|\{1,\psi _n\}|\{1,\psi _n\}\},`$
$`|\{1,n\}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\{|\{1,\psi _n\}+|\{1,\psi _n\}\}.`$
Inserting these equations into (70) and taking some algebras we obtain

$`H_{F}^{}{}_{}{}^{^{\prime \prime }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\stackrel{m,n}{mn}}{}}{\displaystyle \underset{\sigma ,\sigma ^{^{}}}{}}\text{e}^{it\mathrm{\Omega }(mn)}\{m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma }{2}}|\{\sigma ,\psi _m\}\{\sigma ^{^{}},\psi _n\}|+`$ (72)
$`m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma ^{^{}}}{2}}|\{\sigma ,\psi _m\}\{\sigma ^{^{}},\psi _n\}|\}.`$
For simplicity in (71) we set in the following

$$E_{n,\sigma }=\frac{\mathrm{\Delta }}{2}\sigma n|\text{e}^{x(L_+L_{})}|n,$$
(73)
then

$$E_{n,\sigma }=\{\begin{array}{cc}(N)\frac{\mathrm{\Delta }}{2}\sigma \text{e}^{\frac{2g^2}{\omega ^2}}L_n\left(\frac{4g^2}{\omega ^2}\right)\hfill & \\ (K)\frac{\mathrm{\Delta }}{2}\sigma \frac{n!}{(2K)_n}(1+|\kappa |^2)^{Kn}F_n(|\kappa |^2:2K)\text{where}\kappa =\text{sinh}(x)\hfill & \\ (J)\frac{\mathrm{\Delta }}{2}\sigma \frac{n!}{{}_{2J}{}^{}P_{n}^{}}(1|\kappa |^2)^{Jn}F_n(|\kappa |^2:2J)\text{where}\kappa =\text{sin}(x)\hfill & \end{array}$$
(74)
from (59) and the results in sectin 3.1. Now let us solve (66)

$$i\frac{d}{dt}\mathrm{\Psi }_0=\frac{\mathrm{\Delta }}{2}H_F\mathrm{\Psi }_0=\frac{\mathrm{\Delta }}{2}(H_{F}^{}{}_{}{}^{^{}}+H_{F}^{}{}_{}{}^{^{\prime \prime }})\mathrm{\Psi }_0.$$

For that if we set $`\mathrm{\Psi }_0(t)`$ as

$$\mathrm{\Psi }_0(t)=\underset{\sigma }{}\underset{n}{}\text{e}^{itE_{n,\sigma }}a_{n,\sigma }(t)|\{\sigma ,\psi _n\},$$
(75)
then we have a set of complicated equations with respect to $`\{a_{n,\sigma }\}`$, see . But it is almost impossible to solve them. Therefore we make a daring assumption : for $`m<n`$

$$\mathrm{\Psi }_0(t)=\underset{\sigma }{}\text{e}^{itE_{m,\sigma }}a_{m,\sigma }(t)|\{\sigma ,\psi _m\}+\underset{\sigma }{}\text{e}^{itE_{n,\sigma }}a_{n,\sigma }(t)|\{\sigma ,\psi _n\}.$$
(76)
That is, we consider only two terms with respect to $`\{n|n0\}`$. After some algebras we obtain

$`i{\displaystyle \frac{d}{dt}}a_{m,\sigma }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\sigma ^{^{}}}{}}\text{e}^{it(E_{n,\sigma ^{^{}}}E_{m,\sigma })}\text{e}^{it\mathrm{\Omega }(mn)}\left\{m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma }{2}}+m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma ^{^{}}}{2}}\right\}a_{n,\sigma ^{^{}}},`$
$`i{\displaystyle \frac{d}{dt}}a_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\sigma ^{^{}}}{}}\text{e}^{it(E_{m,\sigma ^{^{}}}E_{n,\sigma })}\text{e}^{it\mathrm{\Omega }(nm)}\left\{n|\text{e}^{x(L_+L_{})}|m{\displaystyle \frac{\sigma }{2}}+n|\text{e}^{x(L_+L_{})}|m{\displaystyle \frac{\sigma ^{^{}}}{2}}\right\}a_{m,\sigma ^{^{}}}.`$
But we cannot still solve the above equations exactly (see Appendix), so let us make so–called rotating–wave approximation. The resonance condition is

$$(E_{n,\sigma ^{^{}}}E_{m,\sigma })+(mn)\mathrm{\Omega }=0E_{n,\sigma ^{^{}}}E_{m,\sigma }=(mn)\mathrm{\Omega }$$
(78)
for some $`\sigma `$ and $`\sigma ^{^{}}`$, and we reject the remaining term in (4). Then we obtain simple equations :

Interband Transition Case ($`\sigma \sigma ^{^{}}`$) $`E_{n,\sigma }E_{m,\sigma }=(mn)\mathrm{\Omega }`$

$`i{\displaystyle \frac{d}{dt}}a_{m,\sigma }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}\left\{m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma }{2}}m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma }{2}}\right\}a_{n,\sigma }`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma m|\text{sinh}\left(x(L_+L_{})\right)|na_{n,\sigma },`$
$`i{\displaystyle \frac{d}{dt}}a_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\sigma ^{^{}}}{}}\left\{n|\text{e}^{x(L_+L_{})}|m{\displaystyle \frac{\sigma }{2}}+n|\text{e}^{x(a^{}a)}|m{\displaystyle \frac{\sigma }{2}}\right\}a_{m,\sigma }`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma n|\text{sinh}\left(x(L_+L_{})\right)|ma_{m,\sigma }.`$
Intraband Transition Case ($`\sigma =\sigma ^{^{}}`$) $`E_{n,\sigma }E_{m,\sigma }=(mn)\mathrm{\Omega }`$

$`i{\displaystyle \frac{d}{dt}}a_{m,\sigma }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}\left\{m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma }{2}}+m|\text{e}^{x(L_+L_{})}|n{\displaystyle \frac{\sigma }{2}}\right\}a_{n,\sigma }`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma m|\text{cosh}\left(x(L_+L_{})\right)|na_{n,\sigma },`$
$`i{\displaystyle \frac{d}{dt}}a_{n,\sigma }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{\sigma ^{^{}}}{}}\left\{n|\text{e}^{x(L_+L_{})}|m{\displaystyle \frac{\sigma }{2}}+n|\text{e}^{x(L_+L_{})}|m{\displaystyle \frac{\sigma }{2}}\right\}a_{m,\sigma }`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{2}}\sigma n|\text{cosh}\left(x(L_+L_{})\right)|ma_{m,\sigma }.`$
For simplicity we set

$$=\mathrm{\Delta }n|\text{sinh}\left(x(L_+L_{})\right)|m,^{^{}}=\mathrm{\Delta }n|\text{cosh}\left(x(L_+L_{})\right)|m,$$
(81)
then

$$\mathrm{\Delta }m|\text{sinh}\left(x(L_+L_{})\right)|n=,\mathrm{\Delta }m|\text{cosh}\left(x(L_+L_{})\right)|n=^{^{}}.$$

These are two Rabi frequencies as shown in the following. It is important that Rabi frequencies in our models are given by matrix elements of coherent and generalized coherent operators !

By making use of the results in section 3 and (4), (4), (4) we have

$`(N)\{\begin{array}{cc}=\frac{\mathrm{\Delta }}{2}\sqrt{\frac{m!}{n!}}\left(\frac{2g}{\omega }\right)^{nm}\text{e}^{\frac{2g^2}{\omega ^2}}L_{m}^{}{}_{}{}^{(nm)}\left(\frac{4g^2}{\omega ^2}\right)\{1(1)^{nm}\}\hfill & \\ ^{^{}}=\frac{\mathrm{\Delta }}{2}\sqrt{\frac{m!}{n!}}\left(\frac{2g}{\omega }\right)^{nm}\text{e}^{\frac{2g^2}{\omega ^2}}L_{m}^{}{}_{}{}^{(nm)}\left(\frac{4g^2}{\omega ^2}\right)\{1+(1)^{nm}\}\hfill & \end{array}`$ (84)
$`(K)\{\begin{array}{cc}=\frac{\mathrm{\Delta }}{2}\sqrt{\frac{n!m!}{(2K)_n(2K)_m}}\kappa ^{nm}(1+|\kappa |^2)^{K\frac{n+m}{2}}F_{m}^{}{}_{}{}^{(nm)}(|\kappa |^2:2K)\{1(1)^{nm}\}\hfill & \\ ^{^{}}=\frac{\mathrm{\Delta }}{2}\sqrt{\frac{n!m!}{(2K)_n(2K)_m}}\kappa ^{nm}(1+|\kappa |^2)^{K\frac{n+m}{2}}F_{m}^{}{}_{}{}^{(nm)}(|\kappa |^2:2K)\{1+(1)^{nm}\}\hfill & \end{array}`$ (87)
$`\text{where}\kappa =\text{sinh}\left(x\right)\text{with}x=\text{tanh}^1\left({\displaystyle \frac{2g}{\omega }}\right)`$
$`(J)\{\begin{array}{cc}=\frac{\mathrm{\Delta }}{2}\sqrt{\frac{n!m!}{{}_{2J}{}^{}P_{n}^{}{}_{2J}{}^{}P_{m}^{}}}\kappa ^{nm}(1|\kappa |^2)^{J\frac{n+m}{2}}F_{m}^{}{}_{}{}^{(nm)}(|\kappa |^2:2J)\{1(1)^{nm}\}\hfill & \\ ^{^{}}=\frac{\mathrm{\Delta }}{2}\sqrt{\frac{n!m!}{{}_{2J}{}^{}P_{n}^{}{}_{2J}{}^{}P_{m}^{}}}\kappa ^{nm}(1|\kappa |^2)^{J\frac{n+m}{2}}F_{m}^{}{}_{}{}^{(nm)}(|\kappa |^2:2J)\{1+(1)^{nm}\}\hfill & \end{array}`$ (90)
$`\text{where}\kappa =\text{sin}\left(x\right)\text{with}x=\text{tan}^1\left({\displaystyle \frac{2g}{\omega }}\right)`$
From these we find a constraint between $`m`$ and $`n`$ :

Interband Case$`nm=2N1`$$``$$`n=m+2N1`$ for $`N𝐍`$,

Intraband Case$`nm=2N`$$``$$`n=m+2N`$ for $`N𝐍`$.

Now let us solve (4) and (4).

$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}a_{m,\sigma }\\ a_{n,\sigma }\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}0& \sigma \frac{}{2}\\ \sigma \frac{}{2}& 0\end{array}\right)\left(\begin{array}{c}a_{m,\sigma }\\ a_{n,\sigma }\end{array}\right),`$ (97)
$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}a_{m,\sigma }\\ a_{n,\sigma }\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}0& \sigma \frac{^{^{}}}{2}\\ \sigma \frac{^{^{}}}{2}& 0\end{array}\right)\left(\begin{array}{c}a_{m,\sigma }\\ a_{n,\sigma }\end{array}\right),`$ (105)
so their solutions are given by

$`\left(\begin{array}{c}a_{m,\sigma }(t)\\ a_{n,\sigma }(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\text{cos}(\frac{}{2}t)& i\sigma \text{sin}(\frac{}{2}t)\\ i\sigma \text{sin}(\frac{}{2}t)& \text{cos}(\frac{}{2}t)\end{array}\right)\left(\begin{array}{c}a_{m,\sigma }(0)\\ a_{n,\sigma }(0)\end{array}\right),`$ (112)
$`\left(\begin{array}{c}a_{m,\sigma }(t)\\ a_{n,\sigma }(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}\text{cos}(\frac{^{^{}}}{2}t)& i\sigma \text{sin}(\frac{^{^{}}}{2}t)\\ i\sigma \text{sin}(\frac{^{^{}}}{2}t)& \text{cos}(\frac{^{^{}}}{2}t)\end{array}\right)\left(\begin{array}{c}a_{m,\sigma }(0)\\ a_{n,\sigma }(0)\end{array}\right).`$ (120)
We have obtained some solutions under the rotating–wave approximation. Now it may be suited to compare our results with a recent experimental finding in , but this is beyond our scope. See .

Let us conclude this section by a comment. Our ansatz (76) to solve the equation is too restrictive. We want to use (75) to solve the equation, but it is very hard at this stage.

Problem Find more dynamic methods !

## 5 Quantum Computation

Let us reconsider the results in the preceding section in the light of Quantum Computation. Remind once more that the following arguments are based on the rotating–wave approximation.

Interband Case $`(\sigma =\mathrm{𝟏})`$

$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}0& 0& 0& \frac{}{2}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \frac{}{2}& 0& 0& 0\end{array}\right)\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right).`$ (133)
The solution is

$`\left(\begin{array}{c}a_{m,1}(t)\\ a_{m,1}(t)\\ a_{n,1}(t)\\ a_{n,1}(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}\text{cos}(\frac{}{2}t)& 0& 0& i\text{sin}(\frac{}{2}t)\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ i\text{sin}(\frac{}{2}t)& 0& 0& \text{cos}(\frac{}{2}t)\end{array}\right)\left(\begin{array}{c}a_{m,1}(0)\\ a_{m,1}(0)\\ a_{n,1}(0)\\ a_{n,1}(0)\end{array}\right).`$ (146)
Interband Case $`(\sigma =\mathrm{𝟏})`$

$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& \frac{}{2}& 0\\ 0& \frac{}{2}& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right).`$ (159)
The solution is

$`\left(\begin{array}{c}a_{m,1}(t)\\ a_{m,1}(t)\\ a_{n,1}(t)\\ a_{n,1}(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \text{cos}(\frac{}{2}t)& i\text{sin}(\frac{}{2}t)& 0\\ 0& i\text{sin}(\frac{}{2}t)& \text{cos}(\frac{}{2}t)& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}a_{m,1}(0)\\ a_{m,1}(0)\\ a_{n,1}(0)\\ a_{n,1}(0)\end{array}\right).`$ (172)
Intraband Case $`(\sigma =\mathrm{𝟏})`$

$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}0& 0& \frac{^{^{}}}{2}& 0\\ 0& 0& 0& 0\\ \frac{^{^{}}}{2}& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right).`$ (185)
The solution is

$`\left(\begin{array}{c}a_{m,1}(t)\\ a_{m,1}(t)\\ a_{n,1}(t)\\ a_{n,1}(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}\text{cos}(\frac{^{^{}}}{2}t)& 0& i\text{sin}(\frac{^{^{}}}{2}t)& 0\\ 0& 1& 0& 0\\ i\text{sin}(\frac{^{^{}}}{2}t)& 0& \text{cos}(\frac{^{^{}}}{2}t)& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}a_{m,1}(0)\\ a_{m,1}(0)\\ a_{n,1}(0)\\ a_{n,1}(0)\end{array}\right).`$ (198)
Intraband Case $`(\sigma =\mathrm{𝟏})`$

$`i{\displaystyle \frac{d}{dt}}\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& \frac{^{^{}}}{2}\\ 0& 0& 0& 0\\ 0& \frac{^{^{}}}{2}& 0& 0\end{array}\right)\left(\begin{array}{c}a_{m,1}\\ a_{m,1}\\ a_{n,1}\\ a_{n,1}\end{array}\right).`$ (211)
The solution is

$`\left(\begin{array}{c}a_{m,1}(t)\\ a_{m,1}(t)\\ a_{n,1}(t)\\ a_{n,1}(t)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \text{cos}(\frac{^{^{}}}{2}t)& 0& i\text{sin}(\frac{^{^{}}}{2}t)\\ 0& 0& 1& 0\\ 0& i\text{sin}(\frac{^{^{}}}{2}t)& 0& \text{cos}(\frac{^{^{}}}{2}t)\end{array}\right)\left(\begin{array}{c}a_{m,1}(0)\\ a_{m,1}(0)\\ a_{n,1}(0)\\ a_{n,1}(0)\end{array}\right).`$ (224)
If we can identify (76) with an element in two–qubit space

$$a_{m,1}(t)|00+a_{m,1}(t)|01+a_{n,1}(t)|10+a_{n,1}(t)|11𝐂^2𝐂^2$$
(225)
where $`𝐂^2=\text{Vect}_𝐂\{|0,|1\}`$, then the solutions (146), (172), (198), (224) are kinds of controlled unitary operations (gates) which play a crucial role in Quantum Computation, see for example . For example, (224) is just one of controlled unitary gates expressed graphically as

We note here that controlled unitary gates above are written down as

$$\text{C–Unitary}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& u_{11}& 0& u_{12}\\ 0& 0& 1& 0\\ 0& u_{21}& 0& u_{22}\end{array}\right)$$
(226)
for

$$U=\left(\begin{array}{cc}u_{11}& u_{12}\\ u_{21}& u_{22}\end{array}\right)U(2).$$

A comment is in order.  (and ) is considering the same subject. However the authors in them treated it in the weak coupling regime, while we treated it in the strong coupling regime. It is very interesting to investigate a deep relation (connection) between them.

## 6 Discussion

One of motivations of this study is to apply our results to Holonomic Quantum Computation developped by Italian group (Pachos, Rasetti and Zanardi) and the author, see , , ,  and  and recent , .

In this theory we usually use the effective Hamiltonian of a single–mode field of Kerr medium

$$H_0=XN(N1),N=a^{}a\text{where}X\text{is a constant}$$
(227)
as a background and the real Hamiltonian is in one–qubit case given by

$$H(z,w)=W(z,w)H_0W^1(z,w)$$
(228)
where $`W`$ is a product of coherent operator $`U(z)`$ and squeezed one $`S(w)`$ in section 2. In the above Hamiltonian $`H_0`$ the zero–eigenvalue is two–fold degenerate whose eigenvectors are $`|0`$ and $`|1`$. We set $`|vac=(|0,|1)`$. Then we can construct a connection form $`𝒜`$ on the parameter space $`\{(z,w)𝐂^2\}`$ as

$$𝒜=vac|W^1dW|vac$$
(229)
from (228) where $`d=dz\frac{}{z}+dw\frac{}{w}`$. By making use of this connection we can construct a holonomy group $`Hol(𝒜)`$ ($`U(2)`$) which is in this case equal to $`U(2)`$. In Holonomic Quantum Computation we use this holonomy group as unitary operations in Quantum Computation. The point at issue is that we use not full property of the Hamiltonian but only property of the zero–eigenvalue.

By the way, the Hamiltonian $`H_F`$ in (67)

$$H_F=\frac{\mathrm{\Delta }}{2}U_{0}^{}{}_{}{}^{1}(\sigma _3\mathrm{𝟏}_L)U_0\text{where}L=(N)\text{or}(K)\text{or}(J)$$

is very similar to (228). This system is always two–fold degenerate. Then a natural question arises :

Problem Is it possible to perform a holonomic quantum computation by combining the systems {(N), (K), (J)} ?

This is a very interesting and challenging problem.

Acknowledgment.
The author wishes to thank Marco Frasca for his helpful comments and suggestions.

Appendix

On Equations (4)

Here let us write down full equations of (4) with matrix equation form :

$$i\frac{d}{dt}\left(\begin{array}{c}𝐚_m\\ 𝐚_n\end{array}\right)=\left(\begin{array}{cc}\mathrm{𝟎}& A\\ A^{}& \mathrm{𝟎}\end{array}\right)\left(\begin{array}{c}𝐚_m\\ 𝐚_n\end{array}\right)$$
(230)
where

$$𝐚_k=\left(\begin{array}{c}a_{k,1}\\ a_{k,1}\end{array}\right)\text{for}k=m,n$$

and

$`AA(t)`$ $`=`$ $`\left(\begin{array}{cc}\frac{^{^{}}}{2}\text{e}^{it\left(E_{n,1}+E_{m,1}+\mathrm{\Omega }(mn)\right)}& \frac{}{2}\text{e}^{it\left(E_{n,1}+E_{m,1}+\mathrm{\Omega }(mn)\right)}\\ \frac{}{2}\text{e}^{it\left(E_{n,1}+E_{m,1}+\mathrm{\Omega }(mn)\right)}& \frac{^{^{}}}{2}\text{e}^{it\left(E_{n,1}+E_{m,1}+\mathrm{\Omega }(mn)\right)}\end{array}\right)`$ (233)
$`=`$ $`\text{e}^{it\mathrm{\Omega }(mn)}\left(\begin{array}{cc}\frac{^{^{}}}{2}\text{e}^{it\left(E_{n,1}+E_{m,1}\right)}& \frac{}{2}\text{e}^{it\left(E_{n,1}+E_{m,1}\right)}\\ \frac{}{2}\text{e}^{it\left(E_{n,1}E_{m,1}\right)}& \frac{^{^{}}}{2}\text{e}^{it\left(E_{n,1}E_{m,1}\right)}\end{array}\right)`$ (236)
because $`E_{k,\sigma }=E_{k,\sigma }`$.

We can give (230) a formal solution by infinite series (called Dyson series in Theoretical Physics). Then we meet secular terms.

For example let us consider the following simple equation :

$$\frac{d}{dt}a=\text{e}^{i\omega t}a\text{with}a(0)=c.$$

The solution is given by

$$a(t)=\{\begin{array}{cc}c\text{exp}\left(\frac{\text{e}^{i\omega t}1}{i\omega }\right)\omega 0\hfill & \\ c\text{e}^t\omega =0\hfill & \end{array}$$

That is, we meet the secular term.

By the way, we have known how to handle (simple) secular terms called Renormalization Group Method (Approach), see  for a general introduction.

Frasca in  has applied this method to the above equation. The conclusion is interesting, but seems to be rather involved. We are now reconsidering his approach. Therefore let us present

Problem Solve this matrix equation completely !