format: text

text:
# FERMILAB-Pub-02/033-E Long Baseline Neutrino Physics: From Fermilab to Kamioka

## 1 Introduction

Measurements of the atmospheric neutrino flux with the Super-Kamiokande detector have shown that the deficit of muon neutrinos depends on the zenith angle as expected for neutrino oscillations . This result combined with many other experimental constraints provides strong evidence for $`\nu _\mu \nu _\tau `$ oscillations, with mixing $`>0.88`$% (90% C.L.) and $`\mathrm{\Delta }M^2`$ in the range $`(15)\times 10^3`$ eV2 .

Efforts to study these oscillations with accelerator-produced neutrinos are in progress, using $`\nu _\mu `$ beams produced from pion decay. The K2K collaboration has detected neutrino interactions in the Super-K detector using a neutrino beam from KEK 250 km distant. They report a deficit of events compared to the expectation with no oscillations . The MINOS experiment is under construction, and will use an Iron/scintillator detector 730 km distant from a neutrino beam produced at Fermilab . These experiments are expected to confirm that the $`\nu _\mu `$ disappearance has the $`L/E`$ dependence expected for oscillations, and improve the precision on $`\mathrm{\Delta }M^2`$.

Measurements of solar neutrino fluxes also provide evidence for neutrino oscillations, usually explained as $`\nu _e\nu _\mu `$ oscillations with $`\mathrm{\Delta }m^2`$ much smaller than the $`\mathrm{\Delta }M^2`$ for atmospheric oscillations. The Large Mixing Angle (LMA) solution is strongly favored by the data, with $`\mathrm{\Delta }m^2`$ in the range $`(325)\times 10^5`$ eV2 . The KamLAND  and Borexino  experiments, currently under construction, are expected to test the LMA solution and better constrain the mixing parameters.

The atmospheric and solar neutrino oscillations can be simultaneously explained with mixing in a three-flavor model: The the three known flavors of neutrinos, $`\nu _e`$, $`\nu _\mu `$, and $`\nu _\tau `$, are mixtures of mass eigenstates $`\nu _1`$, $`\nu _2`$, and $`\nu _3`$. The data require a mass hierarchy, with $`|m_3m_2||m_2m_1|`$, The hierarchy can be “normal”, with $`m_3>m_2`$, or “inverted”, with $`m_2>m_3`$. We identify $`\mathrm{\Delta }M^2=m_3^2m_2^2`$ and $`\mathrm{\Delta }m^2=m_2^2m_1^2`$.

While $`\nu _\mu \nu _e`$ oscillations are mainly associated with solar oscillations, they can also occur at the atmospheric frequency if there is a $`\nu _3`$ component in the $`\nu _e`$ eigenstate. Reactor experiments have placed limits on the size of this mixing, resulting in the limit $`\mathrm{sin}^2(2\theta _{13})<0.1`$ . If $`\mathrm{sin}^2(2\theta _{13})>0`$ can be detected, the phenomenology of $`\nu _\mu \nu _e`$ oscillations becomes very rich:

* Interference between the solar and atmospheric oscillation amplitudes can result in large CP violation .
* Matter effects can modify the oscillation probability, depending on the sign of $`\mathrm{\Delta }M^2`$ .

To study these effects on the oscillation probability, we use the parametrization of Ref. , with the oscillation parameters of Table 1 as defaults. For example, Fig. 1 shows the possible CP asymmetry as a function of $`\mathrm{sin}^2(2\theta _{13})`$, at the peak of the atmospheric oscillations. While the mixing in solar oscillations is greater than in atmospheric oscillations, the solar oscillations develop much more slowly, so the two amplitudes for $`\nu _\mu \nu _e`$ can be comparable at atmospheric baselines. We note that maximal CP asymmetry is a possibility: There could be a large signal for neutrinos and no signal at all for antineutrinos, or vice-versa.

Given the potential for profound insights into fundamental physics, most studies of future neutrino oscillation experiments with conventional $`\nu _\mu `$ beams have focused on $`\nu _\mu \nu _e`$ oscillations. A wide variety of neutrino energies and detector technologies have been considered in Refs. .

The potential of low energy neutrino beams combined with large water-cerenkov detectors has been explored for beams based at JHF , CERN , and Fermilab  and baselines of 100 to 300 km. These low energy approaches have the advantage of very low backgrounds, and have excellent sensitivity to very small $`\nu _e`$ appearance signals. They have the disadvantage that matter effects are too small to disentangle from CP violating effects, but are still large enough that ambiguities from the unknown sign of $`\mathrm{\Delta }M^2`$ affect the ability of the experiment to establish CP violation .

The potential of higher energy neutrino beams with very long baseline ($`L>7000`$ km) has been explored in Refs. . The backgrounds in water cerenkov detectors are larger for higher neutrino energies. However, the $`\nu _\mu \nu _e`$ signal at these baselines can be highly amplified by matter effects, improving the signal-to-background ratio . Also, the amplified signal depends mostly on atmospheric oscillations, with solar oscillations playing a much smaller role than with shorter baselines.

From this summary of short- and long-baseline neutrino oscillation experiments, we see that they can play very complimentary roles. The short-baseline experiment measures a CP-violating combination of atmospheric and solar oscillations. The long-baseline experiment determines the sign of $`\mathrm{\Delta }M^2`$, constrains the matter effects, and constrains the size of the atmospheric $`\nu _\mu \nu _e`$ oscillation alone.

We propose herein an optimized long-baseline experiment to complement a short-baseline experiment. For concreteness, we assume a target detector at the Kamioka site in Japan, where the Super-K detector already exists. The JHF-Kamioka proposal  uses the Super-K detector for the first phase, and for the second phase assumes construction of a new detector with 40 times more fiducial mass. We assume that this detector will be used as the target for two neutrino beams: One from JHF, with a baseline of 295 km, and one from Fermilab, with a baseline of 9300 km. The concept of our proposal does not depend on the details of this choice, and could be adapted for other locations. The main principles of our proposal are:

1. Combine information from very short and very long baselines in order to determine the mass hierarchy and CP-violating oscillation parameters.
2. Match the spectrum of the long-baseline neutrino beam to the energy at which matter effects produce the maximum amplification of the signal.
3. Design the long-baseline neutrino beam to have an energy spectrum with a rapid cut-off above the signal region. Since most backgrounds feed-down from the neutrino energy to a lower visible energy, this reduces the background in the signal region.

## 2 The $`\nu _\mu \nu _e`$ probability at 9300 km

We start in the “leading approximation” which parametrizes oscillations driven by the atmospheric $`\mathrm{\Delta }M^2`$ and neglects oscillations driven by the solar $`\mathrm{\Delta }m^2`$. The $`\nu _\mu \nu _e`$ probability can then be written as

$$P(\nu _\mu \nu _e)=R_m^2\mathrm{sin}^2\theta _{23}\mathrm{sin}^22\theta _{13}\mathrm{sin}^2\left(\frac{1.267\mathrm{\Delta }M^2}{R_m}\frac{L}{E_\nu }\right).$$
(1)
The effects of interaction with matter are included in the parameter $`R_m`$. In vacuum, $`R_m=1`$. For positive $`\mathrm{\Delta }M^2`$, there is a resonant energy at which there is full mixing, with $`R_m=1/(\mathrm{sin}\theta _{23}\mathrm{sin}2\theta _{13})`$. For antineutrinos, this resonance occurs for negative $`\mathrm{\Delta }M^2`$. Since we’re considering the case $`\mathrm{sin}^22\theta _{13}<0.1`$, the oscillation length becomes very long, and for terrestrial baselines we can expand the probability to first order in $`L`$:

$$P(\nu _\mu \nu _e)=\mathrm{sin}^2\theta _{23}\mathrm{sin}^22\theta _{13}\left(1.267\mathrm{\Delta }M^2\frac{L}{E_\nu }\right)^2.$$
(2)
Since the flux is proportional to $`1/L^2`$, we find that the number of expected signal events is independent of distance. Since background rates decline with $`L^2`$, longer baselines will provide better signal to noise. For negative $`\mathrm{\Delta }M^2`$, the signal would be suppressed and unobservable.

Fig. LABEL:fig:fkprob shows $`P(\nu _\mu \nu _e)`$ as a function of neutrino energy for a baseline of 9300 km, $`\mathrm{\Delta }M^2=3.5\times 10^3`$ eV2, and $`\mathrm{sin}^22\theta _{13}=0.01`$. We see that the resonant energy is in the range 8 to 10 GeV. At this energy, the oscillation probability is 20 times the maximum probability occurring in vacuum.

Fig. LABEL:fig:fkprob also compares a numerical calculation  of the exact theory incorporating the density variation of the Earth, with an analytical approximation  using an average density. The agreement is excellent in the energy range of interest, and we use the analytical approximation for the rest of this paper. Refs.  have also shown that using an average density is a good approximation for similar parameters.

Fig. LABEL:fig:pvsdm2 shows $`P(\nu _\mu \nu _e)`$ for different values of $`\mathrm{\Delta }M^2`$. The resonant energy is proportional to $`\mathrm{\Delta }M^2`$, while the maximum oscillation probability is constant.

The analytical approximation of Ref.  can also be used to include non-leading terms. The amplification occurs mainly for the leading terms, so the oscillation probability for a 9300 km baseline is dominated by the leading terms much more than for shorter baselines.

## 3 The neutrino beam

Our starting point for the neutrino beamline design is the NuMI beamline at Fermilab, which is currently under construction. In the NuMI beamline, 120 GeV protons from the Main Injector collide with a target which is designed for an average power of 0.4 MW. Pions produced in the target are focused and charge-selected in a double horn system, and then decay in a 675 m decay tunnel. This produces a beam consisting of mainly $`\nu _\mu `$. Different configurations of the target and horns produce spectra peaking at different energies. There is a small $`\nu _e`$ contamination arising from the decay of kaons produced in the target, and the decay of muons produced in pion decays. The level of this contamination varies from 0.6 to 1%, depending on the configuration.

The accelerator used to inject protons into the Main Injector is the Booster ring, which accelerates protons up to 8 GeV. Proton drivers with power of 1 MW or more are currently being designed . With a proton driver replacing the Booster, the Main Injector would be capable of delivering 2 MW of proton power to target. In the case of the linac option for the proton driver, it may be possible to deliver this power at any proton energy less than 120 GeV .

To aim a neutrino beam from Fermilab to the Kamioka site, the decay tunnel must be angled downwards $`46^{}`$ from the horizontal. At the Fermilab site, the decay tunnel must fit within a 200 m vertical depth in order to be within rock favorable for construction and to avoid a water aquifer. This limits the length of the decay tunnel, and reduces the efficiency for pions to decay and produce neutrinos. The length is further reduced by the vertical space needed to bend the protons down to the target, and the space needed for the target and focusing, and beam dump at the end of the tunnel.

As discussed in Sect. 5, most backgrounds for our measurement feed-down from the neutrino energy to a lower visible energy, while for the signal the entire neutrino energy is visible. Backgrounds in the signal region tend to arise from neutrinos with higher energy. Therefore, our goal for the neutrino beamline is to produce a spectrum which peaks at the energy for the maximum oscillation probability in Fig. LABEL:fig:pvsdm2, and falls off quickly at higher energy. A saw-tooth shape would be ideal. The lack of a high-energy tail might also simplify the radiation shielding requirements relative to NuMI.

The beamline will need to be tuned for the actual value of $`\mathrm{\Delta }M^2`$ when it is better known; for now we assume $`\mathrm{\Delta }M^2=3.5\times 10^3`$ eV2, and aim for a peak energy of 9 GeV.

We have used the NuMI fast beamline simulation to produce a candidate beamline design. For the following, we apply a scaling factor of 1.25 needed to bring the fast simulation into agreement with the full simulation of the medium- and high-energy NuMI configurations . Starting with a standard NuMI configuration, we have made the following changes:

1. Use 60 GeV protons on target instead of 120 GeV. This reduces the cost of the bending magnets, and allows the bend to occur in less vertical space. We assume a field of 6 T for the bending magnets.
2. We assume the accelerator complex is upgraded to provide 2 MW of proton power from the Main Injector.
3. Change the target and horn configuration to produce a spectrum peaking at 10 GeV.
4. Add a small dipole bend after the horns, to bend 20 GeV pions by 0.5. This allows the higher energy pions to be filtered out, and reduces some of the $`\nu _e`$ component by filtering out neutral kaons.

The resulting spectrum is shown in Fig. 4, and has the basic properties we were aiming for. The $`\nu _\mu `$ spectrum including oscillation effects is also shown. In an interesting coincidence, the 2nd disappearance peak happens to occur at the peak of the spectrum. The length of the pion decay tunnel is 210 m, and the decay efficiency is 37% what it would be for a NuMI length tunnel.

The event rates do not depend critically on the specifics of the design. For example, if we need to use 120 GeV protons to get the full beam power, the pion decay efficiency relative to NuMI decreases only to 35%. Or, if we use 2 T bending magnets and 60 GeV protons, the relative efficiency decreases to 33%.

Ref.  has also presented beamline simulations producing energy spectra with sawtooth shapes, using quadrupole magnets instead of horns to focus the pions.

### 3.1 Flux normalization

Since the use of space for a near detector would further reduce the efficiency of the beamline, we have assumed that there will be no near detector and the flux can be normalized using data in the far detector. As shown in Fig. 4, we expect a large rate of $`\nu _\mu `$ charged-current interactions in the far detector, and this can in principle be used to normalize the flux. The systematic limitations of this procedure will have to be studied with a detailed simulation. It’s possible that neutral-current interactions or $`\nu _\tau `$ charged-current interactions might also be used to help measure the flux.

### 3.2 Running with antineutrino beams

The neutrino beamline can be converted into an antineutrino beamline by reversing the polarity of the focusing horns and bending magnet. Since the antineutrino cross-section is smaller than the neutrino cross-section, and $`\pi ^{}`$ production in the proton target is suppressed relative to $`\pi ^+`$ production, event rates in the far detector are typically 1/3 of those in neutrino beams. The intrinsic $`\nu _e`$ and $`\overline{\nu }_e`$ backgrounds tend to be worse as well. These problems are mitigated by the following effects:

* As discussed in Sect. 4, the efficiency of reconstructing $`\overline{\nu }_e`$ events is higher than for $`\nu _e`$ events. After reconstruction and event selection cuts, the antineutrino rate will be about 55% of the neutrino rate.
* The dipole bend helps remove neutral kaons and $`K^+`$ in the pion decay tunnel, removing the $`\nu _e`$ component of the beam. The remaining $`\overline{\nu }_e`$ component will be comparable to the $`\nu _e`$ component of the neutrino beam.

## 4 The $`\nu _e`$ signal

The $`\nu _e`$ signal is detected from charged-current (CC) interactions producing an electron along with a hadronic shower with energy $`E_{\mathrm{had}}`$. The visible energy of the event, $`E_{\mathrm{vis}}`$, (the energy not carried away by neutrinos) is equal to the neutrino energy. The differential cross-section in terms of the variable $`y=E_{\mathrm{had}}/E_\nu `$, normalized to 1, can be approximately parametrized by :

$$\frac{d\sigma _\nu }{dy}=\frac{15}{16}\left(1+\frac{(1y)^2}{5}\right).$$
(3)
For antineutrinos, the parametrization is:

$$\frac{d\sigma _{\overline{\nu }}}{dy}=\frac{15}{8}\left(\frac{1}{5}+(1y)^2\right).$$
(4)
We will model the efficiency of detecting a $`\nu _e`$ CC event by assuming the fiducial mass of the detector for the target mass, and by requiring the production of an electron above some threshold $`E`$. The probability of a $`\nu _e`$ with energy $`E_\nu `$ generating an electron above a threshold $`E`$ is:

$$p_{cc}=\frac{15}{16}\left(\frac{16}{15}\frac{E}{E_\nu }\frac{1}{15}\left(\frac{E}{E\nu }\right)^3\right).$$
(5)
For antineutrinos, the parametrization is:

$$p_{cc}=\frac{15}{8}\left(\frac{8}{15}\frac{1}{5}\frac{E}{E_\nu }\frac{1}{3}\left(\frac{E}{E\nu }\right)^3\right).$$
(6)
For neutrinos, we choose the following cuts for the threshold $`E`$:

1. $`E>4.5`$ GeV.
2. $`E/E_{\mathrm{vis}}>0.45`$.

Fig. LABEL:nue\_spectrum shows the $`\nu _e`$ CC energy spectrum at the Kamioka site, assuming the neutrino beam described in Sect. 3, and the neutrino oscillation parameters of Fig. LABEL:fig:fkprob. The dashed line of Fig. LABEL:nue\_spectrum shows the effect of the above cuts. The average efficiency of these cuts is 50%. With these cuts, we find 0.21 events in the signal region, per kton-year of exposure, where the signal region is defined as energy between 7 and 10.5 GeV.

The $`y`$-distribution for antineutrinos is more peaked towards zero than for neutrinos, as described by Eqs. (3) and (4). Therefore, for the same NC background fraction, $`\overline{\nu }_e`$ CC events can be selected with more efficient cuts. For neutrinos, we choose the following cuts for the threshold $`E`$:

1. $`E>3.3`$ GeV.
2. $`E/E_{\mathrm{vis}}>0.33`$.

Fig. LABEL:nueb\_spectrum shows the $`\overline{\nu }_e`$ CC energy spectrum at the Kamioka site, assuming the neutrino beam described in Sect. 3, and the neutrino oscillation parameters of Fig. LABEL:fig:fkprob with negative $`\mathrm{\Delta }M^2`$. The dashed line of Fig. LABEL:nueb\_spectrum shows the effect of the above cuts. The average efficiency of these cuts is 80%, producing an antineutrino signal rate 55% of the neutrino rate.

## 5 Backgrounds

### 5.1 Backgrounds from neutral-current events

Neutral pions are produced in the hadronic showers present in all neutrino interactions with nuclei. At high energies, the photons from the $`\pi ^0\gamma \gamma `$ decay produce a single electromagnetic shower, which in water cerenkov detectors is indistinguishable from that of an electron. Neutral-current (NC) events with an energetic $`\pi ^0`$ constitute a major background to the $`\nu _e`$ appearance signal. The visible energy of these events is given by $`yE_\nu `$, so they feed-down from the neutrino energy to a lower visible energy.

We have studied the production of neutral pions with the LEPTO program  which generates deep inelastic neutrino interactions and models the fragmentation of the hadronic component. There are several efforts underway to model resonance and quasi-elastic production, which we have not included, and which might account for roughly 15% of the cross-section in the energy range we are interested in . However, we do expect the backgrounds to come more from high-$`y`$ events, or from feed-down from higher-energy neutrinos, where the fragmentation assumption is more applicable.

We find that Eq. (3) and Eq. (4) agree well with LEPTO-generated data for charged-current and neutral-current events.

The $`\pi ^0`$ fragmentation functions are roughly independent of hadronic energy. This is illustrated in Fig. 7, which shows integral plots of the $`E(\pi ^0)/E(\mathrm{had})`$ distribution. The probability for a hadronic shower to produce a $`\pi ^0`$ with energy greater than a fraction $`x`$ of $`E(\mathrm{had})`$ can be parametrized by:

$$p(x)=(0.49)(0.96)x+(0.47)x^2.$$
(7)
For a given neutrino energy, we can use Eq. 3 to obtain the visible energy spectrum in NC interactions. For a given visible energy, we can then use Eq. 7 to calculate the probability that there will be a $`\pi ^0`$ above some cut-off energy. Thus, we can convolute an input neutrino spectrum into an output NC background spectrum.

Fig. LABEL:nue\_and\_nc shows the NC background to the $`\nu _e`$ signal, estimated with this procedure with the cuts of Sect. 4. We find a background of 0.047 events per kton-year in the signal region, corresponding to $`f_B=2.1\%`$, where $`f_B`$ is the background fraction relative to the expected signal after cuts for a beam with 100% $`\nu _e`$ . Fig. LABEL:nueb\_and\_nc shows the same distributions for the $`\overline{\nu }_e`$ signal. We obtain the same background fraction with a signal rate 55% of the $`\nu _e`$ rate.

Since the direction of the incoming neutrino is known, and the outgoing neutrino can carry away transverse energy, the direction of the observed energy flow might also be used to further reduce NC background. A detailed simulation will be needed to investigate this possibility.

### 5.2 Background from $`\nu _\mu `$ charged-current interactions

In this analysis, we assume that all $`\nu _\mu `$ CC events can be rejected due to the presence of a muon, and only consider the $`\pi ^0`$ background from NC events. For high-$`y`$ events however, the muon will have low energy and can easily be missed, while the hadronic shower has most of the energy of the neutrino. Neutral pions from these hadronic showers will generate backgrounds the same as for neutral currents. Fortunately, as illustrated in Fig. 4, the $`\nu _\mu `$ flux is very low in the signal region. This is true for other values of $`\mathrm{\Delta }M^2`$ as well. So, we expect this source of background to be negligible compared to the background from neutral current events.

### 5.3 Background from $`\nu _\tau `$ charged-current interactions

At the peak energy for the $`\nu _e`$ signal, most of the $`\nu _\mu `$ have oscillated into $`\nu _\tau `$, so $`\nu _\tau `$ CC interactions are a potentially serious source of background. However, at this energy the $`\nu _\tau `$ CC cross-section is only $``$ 8% of the $`\nu _e`$ CC cross-section. The decay mode $`\tau e\nu \overline{\nu }`$ produces a real electron background, but the branching ratio is only 18%. Neutral pions in $`\nu _\tau `$ CC events can also produce backgrounds. Since all these events have an escaping neutrino from the $`\tau `$ decay, these backgrounds feed-down to lower visible energy. As long as we keep the high energy tail of the neutrino beam suppressed, we expect this source of background to be small. However, the cuts needed to suppress this background need more study with a detailed simulation.

### 5.4 Intrinsic $`\nu _e`$ backgrounds

As discussed in Sect. 3 the neutrino beam has a small amount of $`\nu _e`$ contamination. The dipole bend in the beamline can help reduce this background by eliminating neutral kaons in the beamline. Also, the $`\nu _e`$ energy spectrum is very broad, while the signal region is narrow. The intrinsic $`\nu _e`$ background as a function of the energy width of the beam has been studied in Ref. , showing that it is greatly reduced for narrow beams. Therefore, we expect this background also to be very small.

## 6 Comparison with the JHF-Kamioka project

The first stage of the JHF-Kamoika project uses a 0.75 MW proton driver to generate a neutrino beam aimed at the Super-Kamiokande detector. In phase 2, the proton power is increased to 4 MW, and a new water cerenkov detector, called Hyper-K, is built with 40 times more fiducial mass.

With very similar mixing parameters as those used in Sect. 4, the JHF-Kamioka proposal obtains 0.11 events per kton-year in phase 1, compared to 0.21 for our proposal. The event rate for the Fermilab beam is 200% of the rate for the phase 1 JHF beam, and 40% of the rate for the phase 2 JHF beam.

The JHF-Kamioka project achieves a very low background rate: $`f_B=0.51\%`$, compared to $`f_B=2.5\%`$ for the Fermilab beam, where $`f_B`$ is the background fraction relative to the expected signal after cuts for a beam with 100% $`\nu _e`$ . However, the signal from the Fermilab beam is amplified by a factor of 20, so the effective background for comparison with the JHF beam is only 0.13%.

The signal region for the Fermilab beam is fairly narrow, between 7 and 10.5 GeV. The signal region for the JHF beam is defined as energy between 0.4 and 1.2 GeV. The average oscillation probability over the JHF spectrum is only 56% of the peak probability. This further accentuates the background advantage of the Fermilab beam.

## 7 Combined analysis of data from Fermilab and JHF beams

In the previous section, we showed that the Fermilab and JHF beams would have comparable performance. The JHF phase 2 beam has an advantage in rate, while the Fermilab beam has an advantage in background. However, the real advantage comes in combining information from the two beams. The first step will be the detection of $`\nu _\mu \nu _e`$ oscillations and determination of the sign of $`\mathrm{\Delta }M^2`$. The second step will be constraints on sub-leading effects such as CP violation.

### 7.1 Detecting $`\nu _\mu \nu _e`$ oscillations and the sign of $`\mathrm{\Delta }M^2`$

Fig. LABEL:fig:t13\_reach\_1 shows our estimates of the reach in $`\mathrm{sin}^2(2\theta _{13})`$ for JHF and Fermilab beams with the Super-K detector, assuming positive $`\mathrm{\Delta }M^2`$, and atmospheric oscillations only. The JHF curves are in reasonable agreement with curves in the JHF-Kamioka proposal . If a signal appears with the JHF beam, it should be confirmed by the Fermilab beam within one year of running, otherwise $`\mathrm{\Delta }M^2`$ must be negative. It should then be possible to find the signal with the Fermilab beam in antineutrino mode. The performance of the Fermilab beam is excellent: 5 years of running yields a 3$`\sigma `$ reach of 0.003 in $`\mathrm{sin}^2(2\theta _{13})`$.

Fig. LABEL:fig:t13\_reach\_2 shows a similar plot for the Hyper-K detector. A 5 year run should produce a 3$`\sigma `$ reach of 0.001 or better, with a strong dependence on the background level and systematic uncertainty.

It is interesting that if $`\mathrm{sin}^2(2\theta _{13})`$ is small enough for Eq. (2) to be valid, the long-baseline oscillation probability does not depend explicitly on the size of the matter effect for neutrinos at the resonant energy. If the peak energy for the matter resonance can be determined, the theory of matter effects can be tested and systematic uncertainties from matter effects can be kept small .

### 7.2 Constraints on the CP violating phase

The JHF-Kamioka phase 2 program assumes 2 years of neutrino running, and 6 years of antineutrino running. They have assumed the background will be very well constrained, and for our analysis we have assumed a 2% systematic uncertainty on the background. Except for the highest allowed values of $`\mathrm{sin}^2(2\theta _{13})`$, it will be necessary to either find a way to further reduce the background, or to determine the background to this level. For the Fermilab beam, we have assumed a 5% systematic uncertainty on the background, and an 8 year run.

We have neglected the systematic uncertainty on the flux prediction, but note that for the highest allowed values of $`\mathrm{sin}^2(2\theta _{13})`$, it will need to be constrained to 1%.

For the JHF beam, we have assumed the signal to background ratio in the antineutrino beam is the same as for the neutrino beam, and the 6 year exposure for the antineutrino beam is equivalent to 2 years with the neutrino beam. Without a dipole bend, it is likely that the backgrounds with the antineutrino beam will in fact be much worse. As an alternative, we have considered running for 8 years with the JHF neutrino beam.

We have assumed the parameter values in Table 1, and $`\mathrm{sin}^2(2\theta _{13})`$ $`=`$ 0.01, 0.05, or 0.001. We have not included the effect of uncertainties in parameters other than $`\mathrm{sin}^2(2\theta _{13})`$ and $`\delta `$.

We have assumed that the JHF beam measures oscillation probabilities at 0.8 GeV, and the Fermilab beam at 9.0 GeV, although in reality a more sophisticated energy-dependent analysis will be performed. The expected oscillation probability measurements are shown in Table 2.

Constraints in the $`\mathrm{sin}^2(2\theta _{13})`$-$`\delta `$ plane from various combinations of these measurements are shown in Figs. LABEL:fig:tvd\_01\_f2 through LABEL:fig:tvd\_001\_f8\_k2\_a2. For each input value of $`\mathrm{sin}^2(2\theta _{13})`$, maximal CP violation can be detected at $`3\sigma `$ or better, demonstrating a sensitivity over an impressive range of parameter space. This is possible since, as illustrated in Fig. 1, the maximum CP asymmetry is larger if $`\mathrm{sin}^2(2\theta _{13})`$ is smaller, so the larger asymmetry provides statistical compensation for the smaller signal. Over most of this range, 2 years with the Fermilab and JHF beams combined is enough to reach the $`3\sigma `$ criterion. The constraints are improved with an 8 year run during which the JHF beam runs in antineutrino mode for 6 years. If the JHF antineutrino beam does not provide adequate performance, an 8 year run with Fermilab and JHF neutrino beams provides constraints almost as effective.

## 8 Summary and conclusions

We have investigated the physics potential of very long baseline experiments designed to measure $`\nu _\mu \nu _e`$ oscillation probabilities. The principles of our design are to tune the beam spectrum to the resonance energy for the matter effect, and to have the spectrum cut off rapidly above this energy. The matter effect amplifies the signal, and the cut-off suppresses backgrounds which feed-down from higher energy. The signal-to-noise ratio is potentially better than for any other conventional $`\nu _\mu `$ beam experiment.

The measurements from long- and short-baseline experiments are very complementary. The short-baseline experiment measures a CP-violating combination of atmospheric and solar oscillations. The relatively small matter effects can still confuse the interpretation of the measurements. At long baselines, solar oscillations play a relatively much smaller role. The long-baseline experiment determines the sign of $`\mathrm{\Delta }M^2`$, constrains the matter effects, and constrains the size of the atmospheric $`\nu _\mu \nu _e`$ oscillation alone.

In particular, we have investigated the capabilities of a neutrino beam from Fermilab aimed at the Super-K detector in Japan, 9300 km distant. At this baseline, neutrinos at the resonance energy have an oscillation probability 20 times the maximum in vacuum. We have used a simple model to estimate the signal efficiency and background, and it will be important to check this with a full simulation. For a normal mass hierarchy, a five year run can detect $`\mathrm{sin}^2(2\theta _{13})`$ at 3$`\sigma `$ for values as low as 0.003. If the mass hierarchy is inverted, the beam can be run in antineutrino mode, with a similar signal-to-noise ratio, and event rate 55% as high as for the neutrino mode.

In an interesting coincidence for the 9300 km baseline, the second peak of maximum $`\nu _\mu `$ disappearance occurs at the same energy as the peak for $`\nu _e`$ appearance. Thus, interesting measurements of the parameters of $`\nu _\mu `$ disappearance may also be possible.

Phase 2 of the JHF-Kamioka proposal assumes the construction of the Hyper-K detector with 40 times the fiducial mass of the Super-K detector. Having this detector serve as the target for beams from both JHF and Fermilab is a powerful way of placing constraints in the $`\mathrm{sin}^2(2\theta _{13})`$-$`\delta `$ plane. We have investigated values of $`\mathrm{sin}^2(2\theta _{13})`$ ranging from 0.001 to 0.05, and find that an 8 year run can establish maximal CP violation at 3$`\sigma `$ or better throughout this range. If the JHF antineutrino beam has poor performance, using Fermilab and JHF neutrino beams alone is almost as effective in constraining $`\delta `$.

## 9 Acknowledgments

We thank the NuMI collaboration for the use of their fast beamline simulation. In particular, we thank Debbie Harris for help with the use of this software.


text:
# Large-scale gas disk around the radio galaxy Coma ABased on observations with the Westerbork Synthesis Radio Telescope (WSRT) and the Very Large Array (VLA).

## 1 Gas in and around radio galaxies

Many early-type galaxies in the nearby universe appear to have experienced a recent merger/accretion event. In many cases, this has involved at least one gas-rich companion/galaxy that has brought a significant amount of gas into the galaxy. Indeed, H I, FIR and CO observations of elliptical galaxies have demonstrated that many of these objects have an active and interesting cold interstellar medium, often qualitatively (and sometime also quantitatively) similar to that observed in spirals (see Knapp 1999 for a review). Gas is a key factor in determining the evolution of early-type galaxies (e.g. Kauffmann 1996). Among early-type galaxies, gas-rich systems may represent an important phase in the evolution that many elliptical galaxies go through and they give important information on the formation and evolution of these systems.

Apart from being important for the overall evolution of early-type galaxies, there is also compelling morphological and kinematical evidence that the activity in powerful radio galaxies is triggered by galaxy mergers and interactions (e.g. Heckman et al. 1986, Smith & Heckman 1989, De Koff et al. 1996, Verdoes Klein et al. 1999, Capetti et al. 2000; see e.g. Wilson 1996 for a review). This is also supported by the theoretical results (Kauffmann & Haehnelt 2000) that suggest that the evolution of supermassive black holes is strongly linked to the hierarchical build-up of galaxies. Although these processes are likely to be more efficient and frequent at high redshifts, they are also observed in relatively “nearby” radio galaxies.

If mergers/accretions are indeed important in the evolution/formation of radio galaxies, the question arises what the connection is between their gas properties and those of other early-type galaxies, and whether we can learn something from such a possible connection about the processes that are involved in the AGN activity (Kauffmann & Haehnelt 2000). Although a general connection merger-AGN activity is reasonably well established, a number of issues is still open, like, for example, 1) is the activity triggered by major mergers between gas-rich galaxies or by minor accretion events, 2) what is the relationship with other types of merging systems (e.g. the ultra-luminous infra-red galaxies), 3) at what stage of the merger do the jets and the associated activity occur, and 4) are all early-type galaxies likely to go through phases of radio activity?

Powerful radio galaxies are also frequently associated with kpc-sized optical emission line nebulosities, extending up to tens of kpc from the nucleus. In particular, the ionization of the gas is believed to be due to either anisotropic UV radiation from the active nucleus or to shocks induced by the radio jet and/or the radio lobes. It is likely that in many cases both mechanisms need to be invoked to explain the emission-line properties. One important element that may help to distinguish between the relative importance of the two mechanisms would be to know the distribution of the neutral gas in and around the radio galaxies rather than only the ionization patterns induced by jet or AGN.

It is clear from the above that neutral hydrogen studies can provide a good tool to investigate these different aspects. Usually the H I is too weak to be detected in emission in relatively far away radio galaxies, however certain types of radio galaxies seem to be more promising than others for a possible detection. Among the good candidates is the radio galaxy Coma A. There are several observations that indicate that a complex interaction between the radio structure and a rich interstellar (ISM) medium is occuring in this radio galaxy (Sec. 2). To obtain a more complete knowledge of the ISM in Coma A, we have observed this object in the 21-cm line of HI to try to detect, either in emission or absorption, the gaseous environment of Coma A.

## 2 What is known about Coma A

Coma A (3C277.3) is a well-known radio galaxy ($`z=0.08579`$, Clark 1996) studied in detail by e.g. van Breugel et al. (1985). The radio structure and the radio power ($`\mathrm{log}\mathrm{P}=25.82`$ W Hz-1 at 1.4 GHz) is characteristic of a radio galaxy intermediate between Fanaroff-Riley type I (FRI) and II (FRII), although Coma A is often classified as FRII. In Fig. 1 the radio contours from $`1`$ arcsec resolution VLA observations are shown (see van Breugel et al. 1985 for more details). The radio emission is mainly formed by two “fat” lobes. A bright hotspot is observed only in the northern lobe. The total size is $`100`$ kpc111Throughout this paper we use a Hubble constant of $`H_{}=50`$ km s-1 Mpc-1, and $`q_{}=0`$. At the redshift of Coma A 1 arcsec corresponds to 2.2 kpc.

HST observations (Capetti et al. 2000) show a prominent unresolved optical core and patches of dust. A jet-cloud interaction is occurring in Coma A as shown by the cloud of ionised gas observed next to a radio blob just south of the nucleus. This is clearly seen in HST images (Martel et al. 1999, Capetti et al. 2000, Tadhunter et al. in prep). At this location, the southern radio jet appears to be deflected, it “lights up” and finally it expands into the southern lobe (van Breugel et al. 1985).

A spectacular system of interlocking emission-line arcs and filaments of the same extent as the radio emission has been observed in H$`\alpha `$ (Tadhunter et al. 2000). There is a striking match between the radio emission and the morphology of the ionised gas, in particular the optical filament that curves around the northern radio lobe. This suggests that in Coma A the direct interaction with the radio jet ionises the gas, but that also the expansion of the lobe into the environment plays a role, and vice versa that the gaseous environment determines to some extent the morphology of the radio continuum emission.

## 3 H I absorption against the radio lobe

### 3.1 WSRT observations

We have observed Coma A at the frequency of the redshifted H I line (1308.3 MHz) with the WSRT equipped with the DZB backend using a 10 MHz bandwidth, 128 channels, 2 polarisations and all 14 telescopes. The observations were done in 5 observing sessions, four of them of 12 h (14 Mar 2000, 20 March 2000, 19 Jul 2000) while from the last observation (03 Sep 2000) only 9 h of good data were obtained due to technical problems. This large amount of time is required in order to reach good sensitivity and investigate the presence of H I in emission around the radio galaxy.

The data were calibrated and reduced using the Miriad package (Sault et al. 1995) following standard reduction procedures, with the only exception that particular care was taken in checking the stability of the bandpasses of the antennae. This allowed us to reach a spectral dynamic range sufficiently high (better than a few times $`10^4`$) to be limited by the thermal noise. The data cube we choose to use was made using a robust-Briggs’ weighting (Briggs 1995) with robustness of zero. The noise per channel is $`0.22`$ mJy/beam and the beam size $`38\times 15`$ arcsec with p.a. $`2^{}`$. The spectral resolution is $`18.7`$ km s-1. A continuum image was made using the line-free channels. This image was made using uniform weighting (equivalent to a robust image with robustness = –2) in order to obtain a somewhat higher spatial resolution ($`22\times 11`$ arcsec in p.a. $`2^{}`$).

### 3.2 The H I absorption

The continuum image of Coma A, as made from the WSRT data, is shown in Fig. 2. Given the relatively low spatial resolution, almost all the details visible in Fig. 1 are smoothed out. The two radio lobes are just resolved, while an extended component is also visible. The data clearly show that H I absorption is detected against both lobes at different velocities. This makes Coma A the first radio galaxy known where such absorption is seen at large distances ($``$ 30 kpc) from the nucleus. In the southern lobe, the peak absorption is 4.9 mJy beam-1 centred at a velocity of $`V_{\mathrm{hel}}=25560`$ km s-1. The width of the absorption is FWHM $`95`$ km s-1. In the northern lobe, the peak absorption is 3.6 mJy beam-1 at a velocity of $`V_{\mathrm{hel}}=25780`$ km s-1 and the width is FWHM $`115`$ km s-1. The systemic velocity of Coma A is $`V_{\mathrm{hel}}=25719`$ km s-1 (Clark 1996). Hence, the absorption is blueshifted with respect to the systemic velocity in the southern lobe and redshifted in the northern lobe.

The depth of an absorption line ($`\mathrm{\Delta }S`$) depends on the optical depth ($`\tau `$), the continuum flux density ($`S`$) and the covering factor $`c_\mathrm{f}`$ as $`\mathrm{\Delta }S=c_\mathrm{f}S(1e^\tau )`$. Assuming for the moment a covering factor $`c_\mathrm{f}=1`$, the peak optical depths are $`\tau =0.35\%`$ and $`\tau =0.27\%`$ for the southern and the northern lobe respectively.

The column densities of the neutral hydrogen follow from $`N_{\mathrm{H}\mathrm{I}}=1.8310^{18}T_{\mathrm{spin}}\tau 𝑑v`$ where $`T_{\mathrm{spin}}`$ is the spin temperature in Kelvin and $`v`$ is the velocity in km s-1. Assuming the canonical $`T_{\mathrm{spin}}=100`$ K, we find column densities of $`7.210^{19}`$ atoms cm-2 for the southern lobe and $`6.110^{19}`$ atoms cm-2 for the northern lobe.

### 3.3 H I emission

No H I was detected in emission in the immediate vicinity of Coma A. The 3-$`\sigma `$ upper limit to H I emission in a single resolution element of the cube is $`4.310^{19}`$ atoms cm-2. Over the width of the absorption lines, the 3-$`\sigma `$ upper limit is $`7`$-$`810^{19}`$ atoms cm-2, more or less the column densities seen in absorption in the WSRT data. Thus, if there was extended H I emission associated with Coma A it would have been detected (the VLA data show that the actual column densities are somewhat higher, see Sec. 4). It is therefore likely that the extent of the neutral hydrogen associated with Coma A is limited to the extent of the radio lobes and of the optical emission line gas.

None of the galaxies seen (in projection) near Coma A are detected in emission. At a larger distance H I is detected at the location of a small galaxy at $`\alpha =12^\mathrm{h}54^\mathrm{m}41.8^\mathrm{s}`$ $`\delta =27^{}30^{}15.2^{\prime \prime }`$ (J2000) with $`m_R=18.8\pm 0.35`$ (Odewahn et al. 1995). We derive a systemic velocity for this galaxy of $`V_{\mathrm{hel}}25900`$ km s-1 and the H I profile has a FWHM $`300`$ km s-1. This galaxy is situated $`15^{}`$ from Coma A, corresponding to almost 2 Mpc. Thus, no interaction is expected between this galaxy and Coma A. The H I mass is substantial: $`M_{\mathrm{H}\mathrm{I}}10^{10}`$ $`M_{}`$.

## 4 The spatial distribution of the H I absorption

### 4.1 VLA observations

To be better able to locate where the H I is occurring, we have observed Coma A also with the VLA in its B configuration. These observations were carried out on 15/16 May 2001. In order to have sufficient velocity resolution and to cover a large enough velocity range, we used a bandwidth of 3.125 MHz and 64 channels for each of the two bands available. The two bands were placed with some overlap so that a total band of about 6 MHz is covered. A velocity resolution of 12 km s-1 is obtained using this setup. The observation resulted in 6 hours of good quality data while the data of four hours of observations could not be used due to technical problems.

Also the VLA data were reduced using the Miriad package, again using standard reduction procedures. Initially, two data cubes (one for each band) were made using a robust-Briggs’ weighting of zero (Briggs 1995). This gives a spatial resolution of $`4.5\times 4.0`$ arcsec (p.a. $`70^{}`$). The noise in a single channel is 0.62 mJy beam-1. At this resolution, H I absorption is only weakly detected, especially in the northern lobe. We therefore have produced data cubes by tapering the visibilities in order to improve the sensitivity for more diffuse absorption. For the southern lobe, we have tapered the data with a Gaussian of FWHM $`=8`$ arcsec, resulting in a beam of $`9.4\times 9.0`$ (p.a. 69) and a rms noise of 0.78 mJy beam-1 in a single channel. A Gaussian taper of FWHM $`=10`$ arcsec was necessary to image the region of the H I absorption in the northern lobe. The resulting beam is $`11.5\times 11.0`$ (p.a. 66) and the rms noise is 0.85 mJy beam-1 per channel.

A continuum image was also obtained using the line-free channels taken from one of the observing bands. The noise of the continuum image made with robust-Briggs’ weighting of zero is 0.37 mJy beam-1.

### 4.2 The location of the H I absorption

In Fig. 3 the continuum image as made from the VLA data is shown. Much more detail is visible compared to the image made from the WSRT data. The two radio lobes are now well resolved, while still much of the extended component that is visible in the WSRT image, but not in the image made with the VLA A array, is detected.

In the full resolution VLA data cubes very little absorption is reliably detected. No absorption is detected against the nucleus nor against the location just south of the nucleus where the jet-cloud interaction is occurring (van Bruegel et al. 1985), with upper limits to the optical depth of 4 and 2.7 % respectively.

In the southern lobe the peak absorption is detected only at the 3-$`\sigma `$ level at 2.9 mJy beam-1 in a single channel. This corresponds to a maximum optical depth of $`\tau =1.7\%`$. In Fig. 4 the integrated H I absorption is shown (red contours) at the full resolution. The absorption is found near the region with the brightest continuum emission, with a suggestion that the absorption occurring just SW of this maximum in the continuum. No absorption is detected above the 3-$`\sigma `$ level in the northern lobe. These results indicate that the H I absorption is due to an extended screen in front of the radio lobes.

In order to improve the sensitivity for extended absorption lower resolution data cubes were made as described in Sec. 4.1. In these lower-resolution cubes absorption is detected against both radio lobes. Figure 5 shows where this absorption is detected in the northern radio lobe (at the resolution of the tapered data). The red contours represent the integrated H I absorption. The absorption is somewhat extended and is not occurring near the continuum hotspot. It is clearly displaced from it, in the direction of the nucleus.

In Fig 6 the contours of the H I absorption detected against both radio lobes (at the resolution of the tapered data) are shown, superimposed on the H$`\alpha `$ image.

The peak absorption in the southern lobe is 6.2 mJy beam-1 which corresponds to an optical depth of $`\tau _{\mathrm{max}}=1`$ %. This is higher than in the WSRT data, but this can be explained by the difference in velocity resolution. The width of the absorption profile is FWHM $`70`$ km s-1. For the northern lobe the peak of the absorption is 4.2 mJy beam$`{}_{}{}^{}1`$ ($`\tau _{\mathrm{max}}=0.9`$ %) and FWHM $`60`$ km s-1. The corresponding column densities are $`1.310^{20}`$ atoms cm-2 for the southern lobe and $`8.710^{19}`$ atoms cm-2 for the northern lobe, assuming $`T_{\mathrm{spin}}=100`$ K. By integrating over the area where the absorption is occurring, the estimated H I mass is found to be about $`M_{\mathrm{H}\mathrm{I}}10^9`$ $`M_{}`$.

When comparing with the WSRT results one sees that the values of the peak absorption are, within the errors, the same, but that the widths of the absorption profiles are smaller in the VLA data (see Fig. 7). The VLA observations recover about 70-75% of the absorption detected in the WSRT data. This, and the larger widths of the WSRT profiles implies that more absorption, at a level too faint to be detected in the VLA observations, is present over a larger region of the lobes, at somewhat different velocities. The absorption therefore covers a significant fraction of both radio lobes.

## 5 Neutral and ionised gas: A possible large-scale gas disk

The data discussed above show that there is H I absorption in front of both radio lobes of Coma A. This makes Coma A quite unique: the absorption is not detected against the nuclear regions of a radio galaxy but it is situated at large distances (tens of kpc) from the centre.

The absorption covers a significant fraction of both radio lobes. The strongest absorption is localised in areas at about 15 arcsec ($`33`$ kpc) from the radio core. We do not detect H I absorption against some of the brighter regions of radio emission, e.g. the hot spot in the northern lobe and the bright knot just south of the nucleus.

It is important to understand the relationship between the neutral hydrogen detected at large distances from the centre and the large structure of ionised gas found by Tadhunter et al. (2000). For this it is interesting to compare the velocities of the H I with those of the ionised gas (from Clark 1996). These optical data consist of long-slit spectra taken along the major axis of the ionised-gas distribution. These data do not cover exactly the same regions where the H I absorption is detected, but the difference in location is not too large and a comparison is still meaningful.

The match between the optical velocities and those of the H I is very good (Fig. 8). Both the ionised gas and the neutral gas are blueshifted, by about the same amount, south of the centre and redshifted north of the centre. This close correspondence in velocity between the ionised and neutral gas demonstrates that they are two phases of the same gas structure.

From the morphology of the ionised gas, Tadhunter et al. (2000) argued that Coma A could have experienced a recent merger and the filamentary structure of the ionised gas would represent pre-existing gas structures. The large amount of neutral gas detected in this galaxy via H I absorption, supports the idea that Coma A is likely to have experienced a merger. Considering that a large fraction of the total gas mass must be ionised, this merger must have involved at least several times $`10^9`$ $`M_{}`$ of gas. The merger in Coma A must have involved at least a reasonably sized gas-rich galaxy.

A number of faint galaxies are found within 100 kpc of the nucleus of Coma A and could be associated with the extended H$`\alpha `$ filamentary structures. One of them (see Fig.1 in Tadhunter et al. 2000) is situated (in projection) on the western edge of the southern H I cloud. This galaxy, even if spectroscopically confirmed to be a real companion, appears too small to be the donor of the observed gas because of the large amount of gas detected around Coma A and the large extent over which this gas is found.

One of the aims of the observations presented here was to look for very large-scale H I tails around Coma A. These structures are often observed in recent mergers (e.g. Yun et al. 1994; Hibbard & van Gorkom 1996; Hibbard & Mihos 1995). No such structures have been detected in Coma A. The gas acquired is apparently reasonably settled in an area of a few tens of kpc, instead of the hundred kpc long structures often seen in on-going mergers. Moreover, the kinematics and the more or less symmetric distribution of the gas, suggest that the gas has settled at least to some extent. Considering that this settling, given the size and the velocities observed, would take roughly a few times $`10^8`$ yrs, the merger in Coma A is likely to have happened several times $`10^8`$ yrs ago.

An important constraint to consider in the interpretation is that H I absorption is detected against both radio lobes. Neutral gas must be, therefore, present in front of both of them. Thus, to explain the observed characteristics of both the ionised gas and the neutral hydrogen absorption in Coma A, we have to assume that if the gas is distributed in a disk-like structure, the radio plasma is actually expanding in this disk. The distribution of the ionised gas supports this. The morphology of the ionised gas appears to be distorted by the interaction with the radio lobes, e.g. with shocks fronts sweeping up material into the filamentary structures observed. The fact that the neutral hydrogen seems to show quite broad absorption profiles could be an indication of a kinematical disturbance as the effect of interaction between the radio and the ISM. A first-order anti-correlation is observed between the ionised and the neutral gas, in particular in the northern lobe.

It is also interesting to note that we find the H I in the region where we would expect the gas to be ionised if there is an illuminating quasar. Therefore, the H I measurements provide further evidence for the idea that the quasar illumination in Coma A is weak or absent (Tadhunter et al. 2000).

In many cases where a gas/dust disk is observed in a radio galaxy, the radio axis is perpendicular to this disk. However, this is not a golden rule and applies mainly to small dust lanes near the nucleus (e.g. de Koff et al. 2000, de Ruiter et al. 2001). Moreover, several examples exist, perhaps related to the tri-axial nature of the mass distribution, where the inner gas disk is perpendicular to the outer gas disk (e.g. NGC 5266, Morganti et al. 1997; NGC 2685, Oosterloo et al. (in prep.) and possibly also Centaurus A, Schiminovich et al. 1994). In some of these galaxies the radio jet is directed into the plane of the outer gas disk and interacting with this gas disk. A spectacular example of this is the Seyfert galaxy IC 5063 (Morganti et al. 1998). Hence, a geometry with the lobes expanding into the outer gas disk can be expected and is actually observed in some cases.

Note that in the radio continuum, apart from the two lobes, there are extensions to the east and the west with the western extension larger than the eastern one. The morphology of the ionised gas suggests that the interaction on the eastern side is stronger than on the western side. The larger western extension of the radio continuum on suggests that the gas disk plays a role in confining and shaping the radio continuum structure.

## 6 How unique is Coma A?

Large disks of neutral gas have already been observed in a number of early-type galaxies (see e.g. Oosterloo et al. 2001a,b and ref. therein; van Gorkom 1997 and ref. therein) and the detection of such a disk in Coma A fits naturally in the scenario that these large disks and the AGN activity are merger related. The main difference of Coma A with the galaxies mentioned above is the higher radio continuum power. Such extended gas disks are very difficult to detect if the galaxy is at a relatively large redshift (as most of the radio galaxies are), but our work on Coma A shows that such disks can be studied through H I absorption, an observational technique that can, with current radio telescopes, be used out to redshift $``$3.5.

The only other objects known to us that appears to show characteristics of the neutral gas distribution similar to Coma A are 3C 433 and 3C 234. Recent VLA observations (Morganti et al. 2001) show that at least part of the H I absorption (detected originally with the Arecibo telescope, Mirabel 1989) is situated against the southern radio lobe at about 40 kpc from the nucleus. Unlike Coma A, no ionised gas has been detected in 3C 433 at such a distance from the centre. In 3C 234 H I absorption has been detected toward the NE lobe at about 50 arcsec from the centre (Pihlström, 2001).

On a smaller scale, situations similar to Coma A can be found in a handful of objects. We want to stress that in these cases the scale of the phenomenon is at most a few kpc while in Coma A the absorption is observed at tens of kpc from the nucleus.

One of these objects is PKS 2322–123, the central radio galaxy in the cooling flow cluster Abell 2597 studied by O’Dea et al. (1994) and Taylor et al. (1999). Spatially extended 21-cm H I absorption is detected against this source (O’Dea et al. 1994) indicating that the radio lobes are surrounded by a collection of clouds containing both neutral and ionised gas. Interestingly this object also shows evidence of optical arcs associated with the radio lobes, as observed in the case of Coma A. Koekemoer et al. (1999) suggest that the properties of this source are explained in terms of accretion of gas by the cD during the interaction with a gas-rich galaxy. The same could apply to Coma A, with the main difference that in PKS 2322–123 it occurs at only a few kpc from the nucleus.

Two more examples where the H I absorption is observed in coincidence with the radio lobes instead of the radio core are IC 5063 (Morganti et al. 1998, Oosterloo et al. 2000) and 3C 236 (Conway & Schilizzi 2000). Again, in both cases the distance from the nucleus of the H I absorption is only a few kpc. These are examples of a jet-cloud interaction between the radio plasma and the ISM. It is interesting to note that in IC 5063 (and likely also in 3C 236) the radio axis is parallel to the axis of the large-scale gas disk, as we claim could be the case in Coma A.

More radio galaxies with very extended neutral gas distributions are likely to exist. Extended H I absorption (observed against the Ly$`\alpha `$ emission) has been found in a high fraction of high-$`z`$ radio galaxies (van Ojik et al. 1997). This is considered an indication that high-$`z`$ radio galaxies are located in dense environments and is a diagnostic for the effects of radio jet propagation in this dense medium. Although this phenomenon may be occurring more frequently and more efficiently at high redshifts, in the low redshift radio galaxies described here we may witnessing a similar situation. As the sensitivity of the radio telescopes improves, such extended gas distributions will be found in more radio galaxies. We have recently detected a large disk-like H I structure in the compact radio galaxy B2 0648+27 (Morganti et al. 2001). In this object the radio source is unresolved with arcsecond resolution and therefore H I absorption is only detected against the central source. However, H I in emission is detected out to a radius of $`80`$ kpc.

## 7 Summary

We have detected H I absorption against the radio lobes of the radio galaxy Coma A. This absorption covers a significant fraction of these lobes and occurs at large distances from the nucleus ($`33`$ kpc). The kinematics of the H I matches that of the large structure of ionised gas present in Coma A and the ionised and the neutral gas are two phases of the same structure. No H I is detected in emission outside the area covered by the lobes. The fact that we detect absorption on both sides from the nucleus suggests that the radio lobes are expanding into a large gas disk. The morphology of the ionised gas strongly supports this.

The gas disk appears relatively settled indicating that the AGN activity, assuming the period of such activity lasts at most a few times $`10^7`$ yrs, started at least $`10^8`$ yrs after the merger. The mass of the gas involved in this merger is at least $`10^9`$ $`M_{}`$, suggesting that at least one reasonably-sized gas-rich galaxy was involved.

###### Acknowledgements.

The WSRT is operated by the ASTRON (Netherlands Foundation for Research in Astronomy) with support from the Netherlands Foundation for Scientific Research (NWO). The VLA is a facility of the National Radio Astronomy Observatory, which is operated by Associated Universities Inc., under cooperative agreement with the National Science Foundation.


text:
# New Class of Solvable and Integrable Many-Body Models

## I Introduction

Shortly after the development of the BCS theory of superconductivity1 , the Hamiltonian and the associated computational techniques were adapted to create the pairing model of nuclear structure. This model has been successful in qualitatively accounting for nuclear properties such as the energy gap in the single-particle spectrum of even-even nuclei, odd-even mass differences, and the moments of inertia of deformed nuclei in spite of the fact that the methods used by BCS do not apply to finite systems since they violate number conservation. A great deal of effort was invested into attempts to restore number conservation to the theory. At about the same time, an exact analytical expression2 ; 3  for the eigenstates and eigenvalues of the pairing Hamiltonian with constant interaction matrix elements was published and largely ignored by the nuclear physics community. Shortly thereafter the corresponding many-boson problem was also solved4 . More recently, the attention of the condensed matter community has been directed towards understanding pairing correlations in finite systems by a series of experiments on metallic grains with a size of a few nanometers5 ; 6  and the exact solution has been extensively utilized in studies of these systems6 . A separate and important recent development has been the demonstration that the pairing Hamiltonian, with constant interaction matrix elements, is an integrable system.7  This system is then both solvable and integrable. Another important result is the demonstration of a class of solvable and integrable models involving both pairing forces and an interaction term quadratic in the particle number operators8 ; 9 . This paper draws much of its inspiration from these works and extends their results.

The structure of this paper is as follows. We start in Sec. 2 with a definition of the class of Hamiltonians that we will consider. They are the most general, seniority-conserving Hamiltonians with two-body forces, which describe systems of fermions or bosons. These Hamiltonians have a diagonal single-particle kinetic energy plus a pairing term and a term quadratic in the particle number operators with arbitrary matrix elements in all terms. We then consider a set of operators that have the same structure as the Hamiltonian and derive a set of conditions on the matrix elements in the members of the set that must be satisfied if these operators are required to commute with the Hamiltonian. We obtain an explicit solution for the matrix elements in the two-body part of these constants of the motion in terms of the matrix elements in the Hamiltonian and a set of conditions which the Hamiltonian matrix elements must satisfy. Solutions of these conditions yield integrable Hamiltonians. We show that, due to a symmetry of the conditions, each particular solution gives rise to a large family of solutions. The number of free parameters in each family is greater than the number of single-particle energy levels. We also show that previously published8 ; 9  solvable and integrable models satisfy our conditions. By generating a new particular solution of the conditions, we generate a whole new family of integrable models. In Sec. 3, we show that members of this new family of integrable models are also solvable by giving analytical forms for all their eigenstates and eigenvalues. In Sec. 4, we calculate the occupations of the single-particle levels and some two-body correlation functions in these eigenstates.

## II Hamiltonian and Constants of the Motion

We consider the class of seniority conserving Hamiltonians for Fermi or Bose systems that can be written in the following form $`H=T+U+V`$ , where

$`T`$ $`={\displaystyle \underset{i}{}}2ϵ_in_i`$ (2.1)
$`U`$ $`={\displaystyle \frac{sg}{4}}{\displaystyle \underset{i,j}{}}u_{i,j}\left(\mathrm{\Omega }_i+2sn_i\right)\left(\mathrm{\Omega }_j+2sn_j\right)`$
$`V`$ $`={\displaystyle \frac{g}{2}}{\displaystyle \underset{i,j}{}}v_{i,j}\left(b_i^{}b_j+b_j^{}b_i\right)`$
In these expressions, $`i`$ labels a single-particle energy-level of energy $`ϵ_i`$ and $`g`$ is an interaction strength. The interaction matrices $`u`$ and $`v`$ are real and symmetric. The statistics factor $`s=1\left(+1\right)`$ for Fermi (Bose) statistics. The operators $`n`$ and $`b`$ are defined by

Fermions:

The degeneracy of level $`i`$ is $`2\mathrm{\Omega }_i`$ with the factor $`2`$ coming from time-reversal degeneracy and $`\alpha =\pm 1,\mathrm{},\pm \mathrm{\Omega }_i`$ with

$`n_i`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1}{\overset{\mathrm{\Omega }_i}{}}}\left(a_{i,\alpha }^{}a_{i,\alpha }+a_{i,\alpha }^{}a_{i,\alpha }\right)`$ (2.2)
$`b_i`$ $`={\displaystyle \underset{\alpha =1}{\overset{\mathrm{\Omega }_i}{}}}a_{i,\alpha }a_{i,\alpha }`$
Bosons:

There are 2 possibilities for bosons. In the first case, the single-particle states are not time-reversal eigenstates and the above definition applies. In the second case, the single-particle states are time-reversal eigenstates as in the $`𝐤=\mathrm{𝟎}`$ state or the bound states of a potential. We then define the operators

$`n_i`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1}{\overset{\mathrm{\Omega }_i}{}}}a_{i,\alpha }^{}a_{i,\alpha }`$ (2.3)
$`b_i`$ $`={\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \underset{\alpha =1}{\overset{\mathrm{\Omega }_i}{}}}a_{i,\alpha }a_{i,\alpha }`$
The single-particle operators $`a`$ satisfy (Fermi)Bose (anti)commutation relations which can be used to derive the commutation rules for the $`n`$ and $`b`$ operators

$`[n_i,b_j^{}]`$ $`=\delta _{i,j}b_j^{}`$ (2.4)
$`[b_i,b_j^{}]`$ $`=\delta _{i,j}\left(\mathrm{\Omega }_j+2sn_j\right)`$
The effects of degeneracy are expressed explicitly in the presence of the $`\mathrm{\Omega }_j`$ in the commutation relations and the effects of statistics in the presence of the $`s`$ in $`U`$ and in the factors $`\left(\mathrm{\Omega }_j+2sn_j\right)`$ in the Hamiltonian and the commutation relations. The seniority operators are defined by

$$\nu _i=\underset{\alpha =1}{\overset{\mathrm{\Omega }_i}{}}\left(a_{i,\alpha }^{}a_{i,\alpha }a_{i,\alpha }^{}a_{i,\alpha }\right)$$

and they commute with the Hamiltonian. They count the number of unpaired particles in each single particle level. The form of the interaction in (2.1) is chosen so as to mimic that of a Heisenberg XXZ spin chain model after a quasispin substitution for the fermion pair operators.

The eigenstates of such Hamiltonians factor into a product of a state of the unpaired particles times a state of the paired particles. The energy of the state is the sum of the energies of the unpaired particles, which interact only through $`U`$, plus that of the paired particles. The influence of the unpaired particles on the paired particles comes about through through a modification of the degeneracies $`\mathrm{\Omega }_i`$ \- blocking in the case of fermions and antiblocking for bosons. In what follows, we only consider the seniority-zero paired particles and assume that the degeneracies are effective degeneracies which depend upon which state of the system is being considered.

Following Ref. 7, we look for two-body constants of the motion $`K`$ that have a form similar to that of the Hamiltonian, $`K=X+Y+Z`$ , with

$`X`$ $`={\displaystyle \underset{j}{}}2x_jn_j`$ (2.5)
$`Y`$ $`={\displaystyle \frac{sg}{4}}{\displaystyle \underset{j,k}{}}y_{j,k}\left(\mathrm{\Omega }_j+2sn_j\right)\left(\mathrm{\Omega }_k+2sn_k\right)`$
$`Z`$ $`={\displaystyle \frac{g}{2}}{\displaystyle \underset{j,k}{}}z_{j,k}\left(b_j^{}b_k+b_k^{}b_j\right)`$
where $`y`$ and $`z`$ are real and symmetric matrices to be determined in terms of $`x`$ , $`ϵ`$ , $`u`$ , and $`v`$ so that $`K`$ commutes with the Hamiltonian. Since $`x`$ is a set of free parameters, equal in number to the number of degrees of freedom of the system, the solution of the resulting equations will provide a complete set of constants of the motion specified by a choice of $`x`$ and it will be shown that these constants, with different choices for $`x`$ , commute with each other.

The calculation of the commutator $`[H,K]`$ follows from (2.1), (2.4), and (2.5) and some tedious algebra. The result is

$`[H,K]`$ (2.6)
$`=g{\displaystyle \underset{k,l}{}}\left\{\left(\xi _k\xi _l\right)v_{k,l}+\left(e_ke_l\right)z_{k,l}\right\}\left(b_k^{}b_lb_l^{}b_k\right)`$
$`+{\displaystyle \frac{g^2}{2}}{\displaystyle \underset{j,k,l}{}}\left[v_{k,l}\left(y_{j,k}y_{j,l}\right)\left(u_{j,k}u_{j,l}\right)z_{k,l}+\left(v_{j,k}z_{j,l}v_{j,l}z_{j,k}\right)\right]\left[b_k^{}\left(\mathrm{\Omega }_j+2sn_j\right)b_lb_l^{}\left(\mathrm{\Omega }_j+2sn_j\right)b_k\right]`$
where

$`\xi _k`$ $`=x_k{\displaystyle \frac{1}{2}}sgy_{k,k}`$ (2.7)
$`e_k`$ $`=ϵ_k{\displaystyle \frac{1}{2}}sgu_{k,k}`$
Setting the commutator equal to zero yields the equations

$$z_{k,l}=\frac{\left(\xi _k\xi _l\right)}{\left(e_ke_l\right)}v_{k,l},kl$$
(2.8)
and

$$v_{k,l}\left(y_{j,k}y_{j,l}\right)\left(u_{j,k}u_{j,l}\right)z_{k,l}+\left(v_{j,k}z_{j,l}v_{j,l}z_{j,k}\right)=0,kl$$
(2.9)
We assume that $`v0`$ and first consider the case in which $`u=0`$ . Then (2.9), with $`jk`$ , is

$$v_{k,l}\left(y_{k,k}z_{k,k}y_{k,l}\right)+v_{k,k}z_{k,l}=0,kl$$
(2.10)
Without loss of generality, we can require

$$y_{k,k}=z_{k,k}=0$$
(2.11)
Then

$$y_{k,l}=\frac{\left(x_kx_l\right)}{\left(e_ke_l\right)}v_{k,l},kl$$
(2.12)
which requires that $`v_{k,k}`$ be a constant which can be taken to be one. Equation (2.9), with $`j`$ , $`k`$ and $`l`$ distinct, and, using (2.8) and (2.12), becomes

$$\left(v_{k,l}v_{j,k}v_{j,l}\right)\left(\frac{\left(x_jx_k\right)}{\left(e_je_k\right)}\frac{\left(x_jx_l\right)}{\left(e_je_l\right)}\right)=0,jklj$$

which requires $`v_{j,k}=1`$ and is the simple pairing Hamiltonian with constant interaction matrix elements. This is the result of Ref. 7.

We now assume that both $`u0`$ and $`v0`$ . We can solve (2.9), which is antisymmetric in $`k`$ and $`l`$ , for $`y`$ by setting $`j=k`$ which gives

$$v_{k,l}\left(y_{k,k}z_{k,k}\right)v_{k,l}y_{k,l}+\left(u_{k,l}+v_{k,k}u_{k,k}\right)z_{k,l}=0,kl$$
(2.13)
Without loss of generality, we can require

$$y_{k,k}=z_{k,k}=v_{k,k}u_{k,k}=0$$
(2.14)
Then, using (2.8), we have

$$y_{k,l}=\frac{\left(x_kx_l\right)}{\left(e_ke_l\right)}u_{k,l},kl$$
(2.15)
Equation (2.9), with $`j`$ , $`k`$ and $`l`$ distinct, and using (2.8) and (2.15) provides an integrability condition on $`u`$ and $`v`$ . This condition is linear in the $`x^{}s`$ and, equating their coefficients to zero, leads to the single equation

$$\left(xy\right)u(x,z)v(y,z)\left(xz\right)u(x,y)v(y,z)+\left(yz\right)v(x,y)v(x,z)=0,xyzx$$
(2.16)
where we have introduced the notation

$`u_{j,k}`$ $`=u(e_j,e_k)`$ (2.17)
$`v_{j,k}`$ $`=v(e_j,e_k)`$
and write $`(x,y,z)`$ for $`(e_j,e_k,e_l)`$ . Note that this equation is independent of particle statistics and single-particle degeneracies.

Eq. (2.16) has a curious symmetry that allows the generation of a whole class of solutions from any particular solution. Assume that we have a particular solution of (2.16), $`u(x,y)`$ and $`v(x,y)`$ , then there will exist a class of solutions generated from this particular solution of the form

$`u^{}(x,y)`$ $`={\displaystyle \frac{\left(xy\right)}{f\left(x\right)f\left(y\right)}}u[f\left(x\right),f\left(y\right)]`$ (2.18)
$`v^{}(x,y)`$ $`={\displaystyle \frac{\left(xy\right)}{f\left(x\right)f\left(y\right)}}v[f\left(x\right),f\left(y\right)]`$
where $`f`$ is some well behaved function of its argument. This result follows from the fact that the form of Eq. (2.16) is invariant under the transformation

$$uu^{},vv^{},xf$$
(2.19)
Note that iteration of the transformation (2.18) is equivalent to a redefinition of the function $`f`$ so that each class of solutions is closed under this operation. If we require the functions $`f`$ to be monotonic with single-valued inverses, then the transformations have a group structure. In the following we will look for particular solutions of (2.16) which will provide the basis for classes of solutions generated by the transformations (2.18). The members of each class are specified by the set of numbers $`f_i=f\left(e_i\right)`$ so that the number of free parameters in a class is equal to the number of single-particle energy levels plus the number of free parameters in the particular solution.

We first show that recently published8 ; 9  results are solutions of (2.16). The first class of solutions is based upon the simplest solution of (2.16): $`u=v=1`$ . With this solution, $`U`$ is a function of the total number of particles, and the Hamiltonian is equivalent to the simple pairing model Hamiltonian with constant interaction matrix elements. The members of this class, after the transformation (2.18), are the descendents of this Hamiltonian with

$$u(x,y)=v(x,y)=\frac{xy}{f\left(x\right)f\left(y\right)}$$
(2.20)
This class can also be obtained directly from (2.16) by taking $`u=v`$ and dividing (2.16) by $`v(x,y)v(x,z)v(y,z)`$ to get the equation

$$\frac{\left(xy\right)}{v(x,y)}\frac{\left(xz\right)}{v(x,z)}+\frac{\left(yz\right)}{v(y,z)}=0$$
(2.21)
which is linear in $`1/v`$ . The solutions of this equation are $`\left(xy\right)/v(x,y)=f\left(x\right)f\left(y\right)`$ which is (2.20). A second class of solutions is obtained by writing

$$u(x,y)=v(x,y)w(x,y)$$

and obtaining from (2.16)

$$\frac{\left(xy\right)}{v(x,y)}w(x,z)\frac{\left(xz\right)}{v(x,z)}w(x,y)+\frac{\left(yz\right)}{v(y,z)}=0$$
(2.22)
Solutions of this equation are8 ; 9

$`v(x,y)`$ $`={\displaystyle \frac{q\left(xy\right)}{Sinh\left[q\left(xy\right)\right]}}`$ (2.23)
$`w(x,y)`$ $`=Cosh\left[q\left(xy\right)\right]`$
$`u(x,y)`$ $`={\displaystyle \frac{q\left(xy\right)Cosh\left[q\left(xy\right)\right]}{Sinh\left[q\left(xy\right)\right]}}`$
Note that the first class of solutions is just the $`q0`$ limit of the second.

A third and new class of solutions is obtained by writing

$$v(x,y)=\gamma \left(x\right)\gamma \left(y\right)$$
(2.24)
where $`\gamma ^2\left(x\right)`$ is a polynomial in $`x`$ and $`u(x,y)`$ is a symmetric polynomial in $`x`$ and $`y`$ . Substituting this form into (2.16) gives the equation

$$\left(xy\right)u(x,z)\left(xz\right)u(x,y)+\left(yz\right)\gamma ^2\left(x\right)=0$$
(2.25)
which is linear in $`u`$ and $`\gamma ^2`$ . The solution of this equation is

$`u(x,y)`$ $`=1+\gamma _1\left(x+y\right)+\gamma _2xy`$ (2.26)
$`\gamma \left(x\right)`$ $`=\left(1+2\gamma _1x+\gamma _2x^2\right)^{\frac{1}{2}}`$
with $`\gamma _1`$ and $`\gamma _2`$ free parameters. Note that before the transformation (2.19), the first two terms in $`u`$ can be ignored since they result in terms in the Hamiltonian that are either constants, proportional to the total number of particles, or a modification of the single-particle spectrum. However, they are important after the transformation. Note also that this solution reverts back to the first one for the case $`\gamma _2=\gamma _1^2`$ . The case in which there is no constant term in $`u`$ can be handled by making the replacements $`gg\gamma _0`$ , $`\gamma _i\gamma _i/\gamma _0`$ , $`i=1,2`$ , and then taking the limit $`\gamma _00`$ . The net effect of which is the replacement of $`1`$ in (2.26) by $`0`$ .

Further study of Eq. (2.16) has not produced any other simple solutions that are not the result of symmetry transformations applied to the solutions (2.20), (2.23), and (2.26). In the next sections, we consider the eigenstates of the integrable model (2.26) after a symmetry transformation.

## III Eigenstates

It is a semiempirical theorem that the seniority-zero eigenstates of integrable Hamiltonians of the form (2.1), have the form of a product of correlated pairs. This is true for the states of the Hamiltonians with interaction matrix elements given by (2.20) or (2.23)8 ; 9  and we now show, by construction, that it is also true for the solution (2.26).

Assuming that the eigenstates can be taken as a product of correlated pairs, we write, for a typical eigenstate of $`2N`$ particles,

$$|\psi =\underset{r=1}{\overset{N}{}}B_r^{}|0$$
(3.1)
where

$$B_r^{}=\underset{j}{}\varphi _r\left(e_j\right)b_j^{}$$
(3.2)
Then, commuting the Hamiltonian past the product of $`B^{}`$ ’s , we have

$`\left(H2{\displaystyle \underset{p=1}{\overset{N}{}}}E_pE^{\left(0\right)}\right)|\psi `$ $`={\displaystyle \underset{p=1}{\overset{N}{}}}\left({\displaystyle \underset{r=1,\left(rp\right)}{\overset{N}{}}}B_r^{}\right)\left\{[H,B_p^{}]2E_pB_p^{}\right\}|0`$ (3.3)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{p,q=1\left(pq\right)}{\overset{N^{}}{}}}\left({\displaystyle \underset{r=1,\left(rp,q\right)}{\overset{N}{}}}B_r^{}\right)[[H,B_p^{}],B_q^{}]|0`$
where we have written the energy of the vacuum state as

$$E^{\left(0\right)}=\frac{sg}{4}\underset{i,j}{}u_{i,j}\mathrm{\Omega }_i\mathrm{\Omega }_j$$
(3.4)
and where the parameters $`2E_p`$ are pair-energies to be determined. The energy of the state will be given in terms of these parameters. The commutators can be calculated from (2.4) and (3.2) with the results

$$\left\{[H,B_p^{}]2E_pB_p^{}\right\}|0=2\underset{j}{}\left[\left(e_j\frac{g}{2}\underset{k}{}u_{j,k}\mathrm{\Omega }_kE_p\right)\varphi _p\left(e_j\right)+\frac{g}{2}\underset{k}{}v_{j,k}\mathrm{\Omega }_k\varphi _p\left(e_k\right)\right]b_j^{}|0$$
(3.5)
$$[[H,B_p^{}],B_q^{}]=gs\underset{j,k}{}\left\{u_{j,k}\left[\varphi _p\left(e_j\right)\varphi _q\left(e_k\right)+\varphi _p\left(e_k\right)\varphi _q\left(e_j\right)\right]v_{j,k}\left[\varphi _p\left(e_j\right)\varphi _q\left(e_j\right)+\varphi _p\left(e_k\right)\varphi _q\left(e_k\right)\right]\right\}b_j^{}b_k^{}$$
(3.6)
For the new set of integrable models, we take $`u`$ and $`v`$ to be given by (2.26). Then, after the symmetry transformation (2.18),

$`u_{j,k}`$ $`=u(e_j,e_k){\displaystyle \frac{\left(e_je_k\right)}{\left(f_jf_k\right)}}u(f_j,f_k)`$ (3.7)
$`v_{j,k}`$ $`=v(e_j,e_k)=\gamma \left(e_j\right)\gamma \left(e_k\right){\displaystyle \frac{\left(e_je_k\right)}{\left(f_jf_k\right)}}\gamma \left(f_j\right)\gamma \left(f_k\right)`$
$`\varphi _p\left(e_k\right)`$ $`={\displaystyle \frac{\gamma \left(e_k\right)}{e_k\omega _p}}\varphi _p\left(f_k\right)={\displaystyle \frac{\gamma \left(f_k\right)}{f_k\omega _p}}`$
with $`\omega _p`$ a parameter to be determined. For (3.5), we then have

$`\left\{[H,B_p^{}]2E_pB_p^{}\right\}|0`$ $`=2{\displaystyle \underset{j}{}}[(e_j{\displaystyle \frac{g}{2}}{\displaystyle \underset{k}{}}\mathrm{\Omega }_k{\displaystyle \frac{\left(e_je_k\right)}{\left(f_jf_k\right)}}[u(f_j,f_k)\gamma ^2\left(f_k\right)]E_p)`$ (3.8)
$`+{\displaystyle \frac{g}{2}}{\displaystyle \underset{k}{}}\mathrm{\Omega }_k(e_je_k){\displaystyle \frac{\gamma ^2\left(f_k\right)}{\left(f_k\omega _p\right)}}\varphi _p\left(f_j\right)b_j^{}]|0`$
$`=2{\displaystyle \underset{j}{}}\left[e_j+{\displaystyle \frac{g}{2}}{\displaystyle \underset{k}{}}\mathrm{\Omega }_k\left(e_je_k\right){\displaystyle \frac{u(\omega _p,f_k)}{\left(f_k\omega _p\right)}}E_p\right]\varphi _p\left(e_j\right)b_j^{}|0`$
and for (3.6), after algebraic simplification,

$`[[H,B_p^{}],B_q^{}]`$ $`=gs{\displaystyle \frac{u(\omega _p,\omega _q)}{\left(\omega _p\omega _q\right)}}{\displaystyle \underset{j,k}{}}\left(e_je_k\right)\left[\varphi _p\left(f_j\right)\varphi _q\left(f_k\right)\varphi _p\left(f_k\right)\varphi _q\left(f_j\right)\right]b_j^{}b_k^{}`$ (3.9)
$`=2gs{\displaystyle \frac{u(\omega _p,\omega _q)}{\left(\omega _p\omega _q\right)}}{\displaystyle \underset{j}{}}e_j\left[\varphi _p\left(f_j\right)B_q^{}\varphi _q\left(f_j\right)B_p^{}\right]b_j^{}`$
Substituting these results into (3.3) yields

$`\left(H2{\displaystyle \underset{p=1}{\overset{N}{}}}E_pE^{\left(0\right)}\right)|\psi `$ (3.10)
$`=2{\displaystyle \underset{p=1}{\overset{N}{}}}\left({\displaystyle \underset{r=1,\left(rp\right)}{\overset{N}{}}}B_r^{}\right){\displaystyle \underset{j}{}}[e_j(1+{\displaystyle \frac{g}{2}}{\displaystyle \underset{k}{}}\mathrm{\Omega }_k{\displaystyle \frac{u(\omega _p,f_k)}{\left(f_k\omega _p\right)}}gs{\displaystyle \underset{q=1\left(qp\right)}{\overset{N^{}}{}}}{\displaystyle \frac{u(\omega _p,\omega _q)}{\left(\omega _p\omega _q\right)}})`$
$`{\displaystyle \frac{g}{2}}{\displaystyle \underset{k}{}}\mathrm{\Omega }_ke_k{\displaystyle \frac{u(\omega _p,f_k)}{\left(f_k\omega _p\right)}}E_p\varphi _p\left(f_j\right)b_j^{}]|0`$
Setting this equal to zero yields the equations

$$\left[1+\frac{g}{2}\underset{k}{}\mathrm{\Omega }_k\frac{u(\omega _p,f_k)}{\left(f_k\omega _p\right)}+gs\underset{q=1\left(qp\right)}{\overset{N^{}}{}}\frac{u\left(\omega _{p,}\omega _q\right)}{\left(\omega _q\omega _p\right)}\right]=0,p=1,\mathrm{},N$$
(3.11)
$$E_p=\frac{g}{2}\underset{k}{}\mathrm{\Omega }_ke_k\frac{u(\omega _p,f_k)}{\left(f_k\omega _p\right)}$$
(3.12)
The energy of the state is then given by

$`E`$ $`=E^{\left(0\right)}+2{\displaystyle \underset{p=1}{\overset{N}{}}}E_p`$ (3.13)
$`=E^{\left(0\right)}g{\displaystyle \underset{p=1}{\overset{N}{}}}{\displaystyle \underset{k}{}}\mathrm{\Omega }_ke_k{\displaystyle \frac{u(\omega _p,f_k)}{\left(f_k\omega _p\right)}}`$
For $`f_k=e_k`$ , this expression simplifies to

$$E=E^{\left(0\right)}+2\underset{p=1}{\overset{N}{}}\omega _pgs\underset{p,q=1\left(qp\right)}{\overset{N^{}}{}}u\left(\omega _{p,}\omega _q\right)g\underset{p=1}{\overset{N}{}}\underset{k}{}\mathrm{\Omega }_ku(\omega _p,e_k)$$
(3.14)
The equations (3.11) can be transformed into a form that allows a two-dimensional electrostatic analogy for fermions with $`s=1`$ Some algebraic manipulation of (3.11) yields the equations

$$\frac{2+g\left[\mathrm{\Omega }+2s\left(N1\right)\right]\left(\gamma _1+\gamma _2\omega _p\right)}{2g\gamma ^2\left(\omega _p\right)}\frac{1}{2}\underset{k}{}\frac{\mathrm{\Omega }_k}{\left(f_k\omega _p\right)}s\underset{q=1\left(qp\right)}{\overset{N^{}}{}}\frac{1}{\left(\omega _q\omega _p\right)}=0,p=1,\mathrm{},N$$
(3.15)
For fermions, these equations can be interpreted as the equations for the unstable equilibrium positions of $`N`$ unit line charges at the locations $`\omega _p`$ in the complex plane interacting with charges $`\mathrm{\Omega }_k/2`$ at the locations $`f_k`$ and two charges located at the zeros of $`\gamma ^2\left(\omega \right)`$ . In this analogy, the statistical repulsion between fermions is mapped onto the electrostatic repulsion of like charges in two dimensions. Such equations have been shown to lead to the BCS equations plus corrections in a carefully defined $`N\mathrm{}`$ limit.10  Eq. (3.14) then shows that the energy of the states has terms proportional to the dipole and quadrupole moments of the free charge distribution. For bosons, we have a statistical attraction between particles and the electrostatic analogy does not work. Such equations have been shown to yield a generalized Bose condensation into the lowest two single-particle levels for a repulsive interaction in the $`N\mathrm{}`$ limit.4  The Bogoliubov approximation, modified to accommodate this generalized condensation, then gives good results in this limit.

The same techniques can be used to calculate the eigenvalues of the integrals $`K`$ . We consider the solution (2.26) after a symmetry transformation. We then have

$$K|\psi =\left[K^{\left(0\right)}g\underset{p=1}{\overset{N}{}}\underset{k}{}\mathrm{\Omega }_kx_k\frac{u(f_k,\omega _p)}{\left(f_k\omega _p\right)}\right]|\psi $$
(3.16)
where

$$K^{\left(0\right)}=\frac{sg}{4}\underset{j,k}{}\mathrm{\Omega }_j\mathrm{\Omega }_ky_{j,k}$$
(3.17)
If we make the choice $`x_k\delta _{i,k}`$ and $`KK_i`$ , then $`K_i`$ and its eigenvalues, $`k_i`$ , are

$`K_i`$ $`=2n_i+g{\displaystyle \underset{k\left(ki\right)}{}}{\displaystyle \frac{1}{\left(f_if_k\right)}}\left[v(f_i,f_k)\left(b_i^{}b_k+b_k^{}b_i\right){\displaystyle \frac{s}{2}}u(f_i,f_k)\left(\mathrm{\Omega }_i+2sn_i\right)\left(\mathrm{\Omega }_k+2sn_k\right)\right]`$ (3.18)
$`k_i`$ $`=g\left[K_i^{\left(0\right)}\mathrm{\Omega }_i{\displaystyle \underset{p=1}{\overset{N}{}}}{\displaystyle \frac{u(f_i,\omega _p)}{\left(f_i\omega _p\right)}}\right]`$
$`K_i^{\left(0\right)}`$ $`={\displaystyle \frac{s}{2}}\mathrm{\Omega }_i{\displaystyle \underset{k\left(ki\right)}{}}\mathrm{\Omega }_k{\displaystyle \frac{u(f_i,f_k)}{\left(f_if_k\right)}}`$
This result will be used in the next section to calculate the occupations of the single-particle levels and some two-body correlation functions in these eigenstates.

## IV Occupations and Correlations

We can use the eigenvalues of the Hamiltonian and the $`K`$ ,s and the Hellman-Feynman14  theorem to calculate certain two-body correlations in the eigenstates presented in the previous section. The Hamiltonian and constants of the motion are functions of the parameters in the potential, $`g`$ , $`\gamma _1`$ , and $`\gamma _2`$ as well as the single-particle spectrum $`e_j`$ , and its transform $`f_j=f\left(e_j\right)`$ . Since the function $`f`$ is unspecified, we can treat $`e_j`$ and $`f_j`$ as independent parameters. In this picture the $`\omega _p`$ are functions of $`g`$ and the $`f_j`$ , through (3.11), and do not depend upon the $`e_j`$ . Furthermore, comparing (3.13) with (3.18) yields the result

$$k_iK^{\left(0\right)}=\frac{}{e_i}\left(EE^{\left(0\right)}\right)$$
(4.1)
In order to calculate the mean occupations of the single-particle levels, we use (3.18) to write

$$n_i=\frac{1}{2}\left(K_ig\frac{K_i}{g}\right)$$
(4.2)
The Hellman-Feynman theorem then yields

$`\psi \left|n_i\right|\psi `$ $`={\displaystyle \frac{1}{2}}\left(k_ig{\displaystyle \frac{k_i}{g}}\right)`$ (4.3)
$`=\mathrm{\Omega }_i{\displaystyle \underset{p=1}{\overset{N}{}}}{\displaystyle \frac{\gamma ^2\left(f_i\right)}{\left(f_i\omega _p\right)^2}}R_p`$
$`R_p`$ $`={\displaystyle \frac{g^2}{2}}\left({\displaystyle \frac{\omega _p}{g}}\right)`$
The $`R_p`$ satisfy the linear set of equations

$$\left[\underset{k}{}\mathrm{\Omega }_k\frac{\gamma ^2\left(f_k\right)}{\left(f_k\omega _p\right)^2}+2s\underset{q=1\left(qp\right)}{\overset{N^{}}{}}\frac{\gamma ^2\left(\omega _q\right)}{\left(\omega _p\omega _q\right)^2}\right]R_p2s\underset{q=1\left(qp\right)}{\overset{N^{}}{}}\frac{\gamma ^2\left(\omega _p\right)}{\left(\omega _p\omega _q\right)^2}R_q=1$$
(4.4)
obtained by differentiation of (3.11) with respect to $`g`$ .

The Hellman-Feynman theorem can also be used to calculate the expectation values of the following operators

$`{\displaystyle \frac{K_i}{f_j}}`$ $`=g{\displaystyle \frac{1}{\left(f_if_j\right)^2}}{\displaystyle \frac{\gamma \left(f_i\right)}{\gamma \left(f_j\right)}}K_{i,j},ij`$ (4.5)
$`K_{i,j}`$ $`=\left\{u(f_i,f_j)\left(b_i^{}b_j+b_j^{}b_i\right){\displaystyle \frac{s}{2}}\gamma \left(f_i\right)\gamma \left(f_j\right)\left(\mathrm{\Omega }_i+2sn_i\right)\left(\mathrm{\Omega }_j+2sn_j\right)\right\}`$
or

$$K_{i,j}=\frac{\left(f_if_j\right)^2}{2g}\left[\frac{\gamma \left(f_j\right)}{\gamma \left(f_i\right)}\frac{K_i}{f_j}+\frac{\gamma \left(f_i\right)}{\gamma \left(f_j\right)}\frac{K_j}{f_i}\right]$$
(4.6)
The expectation value of $`K_i/f_j`$ is given by

$$\psi \left|\frac{K_i}{f_j}\right|\psi =\frac{k_i}{f_j}$$
(4.7)
and

$`\psi \left|K_{i,j}\right|\psi `$ $`={\displaystyle \frac{\left(f_if_j\right)^2}{2g}}\left[{\displaystyle \frac{\gamma \left(f_j\right)}{\gamma \left(f_i\right)}}{\displaystyle \frac{k_i}{f_j}}+{\displaystyle \frac{\gamma \left(f_i\right)}{\gamma \left(f_j\right)}}{\displaystyle \frac{k_j}{f_i}}\right]`$ (4.8)
$`=\mathrm{\Omega }_i\mathrm{\Omega }_j\gamma \left(f_i\right)\gamma \left(f_j\right)\left[{\displaystyle \frac{s}{2}}+{\displaystyle \frac{\left(f_if_j\right)^2}{4}}{\displaystyle \underset{p=1}{\overset{N}{}}}\gamma ^2\left(\omega _p\right)\left({\displaystyle \frac{R_{p,j}}{\left(f_i\omega _p\right)^2}}+{\displaystyle \frac{R_{p,i}}{\left(f_j\omega _p\right)^2}}\right)\right]`$
where we have written

$$\frac{\omega _p}{f_j}=\mathrm{\Omega }_jR_{p,j}$$
(4.9)
and where the $`R_{p,j}`$ satisfy the linear set of equations (4.4), with the right-hand-side replaced by $`\gamma ^2\left(\omega _p\right)/\left(f_j\omega _p\right)^2`$ , obtained by differentiating (3.11),

## V Conclusions

We have derived a general set of conditions for the integrability of the most general seniority conserving Hamiltonian with two-body forces for bosons or fermions. They have been used to obtain a new class of integrable models. We have shown that the eigenvalues and eigenstates of the Hamiltonian and the associated integrals of the motion can be obtained by the solution of a fairly simple set of coupled algebraic equations. While this solvability depends upon the integrability conditions in a rather complicated way, we have not been able to construct a general result that integrability implies solvability.

The results of this paper can be generalized to other systems such as a system of nucleons with an isospin-independent $`J=0`$ pairing in jj-coupling11  or a spin- and isospin-independent $`L=0`$ pairing in ls-coupling.12  Another system to which these methods can be applied is the system of fermions with an equal mix of pairing in the $`{}_{}{}^{1}S_{0}^{}`$ and $`{}_{}{}^{3}P_{0}^{}`$ channels.13