text:
# Detecting the WIMP-wind via spin-dependent interactions

## 1 Introduction

Existence of the dominant Cold Dark Matter(CDM) became much more concrete by the recent results of the WMAP cosmic microwave background(CMB) all-sky observation with other finer scale CMB measurements (ACBAR and CBI), 2dFGRS measurements, and Lyman $`\alpha `$ forest data . Best fit cosmological parameters, $`h=0.71_{0.03}^{+0.04}`$, $`\mathrm{\Omega }_\mathrm{b}h^2=0.0224\pm 0.0009`$, $`\mathrm{\Omega }_\mathrm{m}h^2=0.135_{0.009}^{+0.008}`$, and $`\mathrm{\Omega }_{\mathrm{tot}}h^2=1.02\pm 0.02`$, show that the CDM consists about 20$`\%`$ of the energy and dominates the mass of the universe.

Weakly Interacting Massive Particles (WIMPs) are one of the best candidates for CDM. In spite of quite a few WIMP search experiments, there has not yet been any experimental evidence of WIMP detection , except for an indication of an annual modulation signal, reported by the DAMA group . Because the amplitude of an annual modulation signal is only a few $`\%`$, positive signatures of WIMPs are very difficult to detect with conventional methods, which basically measure the recoil spectrum. Owing to the motion of the solar system with respect to the galactic halo, the direction-distribution of the WIMP velocity that we observe at the earth is expected to show an asymmetry, like a wind of WIMPs. Attempts to detect a positive signature of WIMPs by measuring the recoil angles have been carried out  ever since it was indicated to be an alternative and more reliable method . Gaseous detectors are one of the most appropriate devices for detecting this WIMP-wind. Properties of the $`\mathrm{CS}_2`$ gas, which is sensitive to the WIMP-wind mainly via spin-independent (SI) interactions are mainly studied because of its small diffusions. We, on the other hand, are focusing on the detection of WIMPs via spin-dependent (SD) interactions, because WIMPs can basically be detected both via SI and SD interactions. Because fluorine was said and also found to be very effective for a WIMP search via SD interactions, as shown in previous works , we studied the detection feasibility with carbon tetrafluoride ($`\mathrm{CF}_4`$) gas, which has been studied very well as one of the standard gases for time projection chambers (TPC) . While the use of $`\mathrm{CF}_4`$ has been proposed by others before , this is the first in-detail feasibility study, taking into account the performance of an existing three-dimensional tracking detector ($`\mu `$-TPC) and the measured neutron background flux.

## 2 $`\mu `$-TPC

A $`\mu `$-TPC is a time projection chamber with a micro pixel chamber ($`\mu `$-PIC) readout, being developed for the detection of tracks of charged particles with fine spatial resolutions . A $`\mu `$-PIC is a gaseous two-dimensional position-sensitive detector manufactured using printed circuit board (PCB) technology. With PCB technology, large-area detectors can, in principle, be mass-produced, which is an inevitable feature for a dark matter detector. We developed a prototype $`\mu `$-TPC with a detection volume of 10 $`\times `$ 10 $`\times `$ 8 $`\mathrm{cm}^3`$, and studied its fundamental properties with a Ar-$`\mathrm{C}_2\mathrm{H}_6`$ gas mixture of 1atm. Anode electrodes are formed with a pitch of 400 $`\mu `$m, and signals from anode strips and cathode strips are amplified and discriminated into LVDS-level signals. Discriminated digital signals are fed to an FPGA-based position encoding module and synchronized with an internal clock of 20 MHz, so that two-dimensional hit positions can be calculated by the anode-cathode coincidence within one clock. Together with the timing information, three-dimensional tracks are thus detected as successive points. Measured two-dimensional ($`\mu `$-PIC intrinsic) and three-dimensional position resolutions of the $`\mu `$-TPC are 120$`\mu `$m and 260$`\mu `$m, respectively. We irradiated the $`\mu `$-TPC with fast neutrons from a radioactive source of $`{}_{}{}^{252}\mathrm{Cf}`$ , and tracks of the recoil protons (500 keV - 1 MeV) were detected, as shown by the filled circles in FIG. 1. The distance between the detected points was restricted by the clock of the electronics (20 MHz), which we will soon increase to more than 50 MHz . This result indicates that a $`\mu `$-TPC, which can detect tracks of the charged particles down to 3 mm, will be available for the dark matter search experiment. The sums of the all cathode signals recorded by a flash ADC (FADC) are also shown in FIG. 1 in the cathode=0 plane. We recognized that Bragg curves were well-reproduced by the FADC data, and thus the track directions were obviously known. In FIG. 1, a typical electron track of about 100 keV is also shown by the blank circles. Because the electron tracks are much more winding, and have a smaller energy deposition ($`dE/dx`$) than those of protons, gamma-rays are known to be discriminated with high efficiency . Stable operation of the $`\mu `$-TPC with a $`dE/dx`$ threshold of 15 keV/cm (gain $``$ 5000) for more than 1000 hours was already realized. We are improving the geometrical structure of the electrodes by using three-dimensional simulators , which we expect will increase the gas gain by a factor of three. Then, the $`dE/dx`$ threshold will be better than 10 keV/cm, even with $`\mathrm{CF}_4`$ gas, which has a larger $`W`$ value (average energy to produce one electron-ion pair) than argon($`W_{\mathrm{CF}_4}`$=54 eV, $`W_{\mathrm{Ar}}`$=26 eV). We started a study of the detector performance with $`\mathrm{CF}_4`$ gas, and found that $`\mu `$-TPC shows a good detector performance (gas gain $``$ 10000). Precise studies on the gas properties will be reported in another paper. Based on the detector studies described above, and on-going improvements, we assume that a $`\mu `$-TPC as a dark matter detector have the properties listed below:

* Track length threshold is 3 mm.
* $`dE/dx`$ threshold is 10 keV/cm.
* Gamma-ray discrimination efficiency is close to 100$`\%`$.
* Bragg curves are detected.

## 3 WIMP-wind measurement

We studied the feasibility of WIMP-wind detection with $`\mathrm{CF}_4`$ gas considering the detector characteristics listed above. First, we calculated the energy deposition of a F ion in $`\mathrm{CF}_4`$ gas. We set the gas pressure at 20 Torr and the recoil energy at 25 keV. Because the binding energy of a C-F in a $`\mathrm{CF}_4`$ molecule is 4.6 eV, all of the recoil energy is thought to become kinetic energy of the F ion. The energy deposition was calculated by SRIM2003  and the range was scaled by the same manner as discussed in Ref. . The energy deposition was also scaled by the same scaling factor. The range $`R`$ of the F ion in $`\mathrm{CF}_4`$ is expressed as a function of the pressure $`p`$ and the recoil energy $`E`$ with a good approximation down to $`E`$=5 keV in this situation as follows:

$$R=\left(\frac{dE}{dx}\right)_{\mathrm{Ne}_{\mathrm{in}}\mathrm{Ar}}^1\frac{N_{\mathrm{Ar}}}{N_{\mathrm{CF}_4}}E\frac{760\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{r}}{p}=kEp^1,$$
(1)
where $`\left(\frac{dE}{dx}\right)_{\mathrm{Ne}_{\mathrm{in}}\mathrm{Ar}}`$ $``$ 130keV/mm  is the energy deposition of a Ne ion in Ar gas; $`\frac{N_{\mathrm{Ar}}}{N_{\mathrm{CF}_4}}`$ is the ratio of the electron density of Ar gas and $`\mathrm{CF}_4`$ gas. We consequently knew that a 25 keV F ion has a range of roughly 3 mm in 20 Torr $`\mathrm{CF}_4`$ using k=2.63 $`\mathrm{mm}\mathrm{Torr}/\mathrm{keV}`$. The calculated energy losses of a F ion in $`\mathrm{CF}_4`$ gas are shown in FIG. 2. The range of an ion is determined by the total (electron + nuclear fields) energy loss, shown in the dashed line. Since the energy loss in the nuclear field has a negligible effect on ionizations, the ”visible” energy loss is that in the electron field shown by the solid line. From the calculated energy loss in the electron field, it is seen that the electron diffusion should be suppressed to lower than $`\sigma <`$ 1 mm in order to know the track direction.

In the following discussion, we consider two experimental modes:

1. We use the digital hit points for tracking, but do not use the FADC data for determining the track direction. In this case, the recoil angle is known only by the absolute value of its cosine. We hereafter call this semi-tracking, or the ”ST” mode.
2. We somehow suppress the electron diffusion and use the full information on their tracks, including the direction. We hereafter call this full-tracking, or the ”FT” mode.

The FT mode will be realized by limiting the drift length, or by using an elencro-negative gas . A magnetic field doesn’t help, however, because the longitudinal electron diffusion cannot be suppressed by the magnetic field. The electron diffusion ($`\sigma `$) along the track path is expressed as $`\sigma =\sqrt{2DL/W}`$, where D is the diffusion coefficient, $`L`$ is the drift length, and $`W`$ is the drift velocity. $`D<`$ 1000 $`\mathrm{cm}^2/\mathrm{s}`$  and $`W>`$ 10 cm/$`\mu `$s are achieved for a reduced electric field of $`E/p>1`$kV/cm/atm, which implies that $`\sigma <`$ 1 mm is realized with a drift length of 50 cm in 20 Torr of $`\mathrm{CF}_4`$ gas. We calculated the recoil angle distributions for the ST and FT modes. We followed Ref.  for the energy-spectrum calculation and

$$\frac{dR}{dEd\mathrm{cos}\gamma }\mathrm{exp}\left[\frac{\left(v_\mathrm{s}\mathrm{cos}\gamma v_{\mathrm{min}}\right)^2}{v_0^2}\right]$$
(2)
in Ref.  for the direction distribution calculation. Here, $`R`$ is the count rate, $`E`$ is the recoil energy, $`\gamma `$ is the recoil angle, $`v_\mathrm{s}`$ is the solar velocity with respect to the galaxy, $`v_{\mathrm{min}}`$ is the minimum velocity of WIMPs that can give a recoil energy of $`E`$, and $`v_0`$ is the Maxwellian WIMP velocity dispersion. We used the astrophysical, nuclear, and experimental parameters given in Table 1. We assumed that the gamma-rays, $`\beta `$-rays, and $`\alpha `$-rays are discriminated by 100$`\%`$ efficiency, and that fast neutrons dominate the background. This assumption is reasonable because we can detect three-dimensional tracks with fine samplings in contrast to the two-dimensional tracks discussed in Ref. . The measured fast neutron flux at Kamioka Observatory  was used for the calculation.

The simulated $`\mathrm{cos}\gamma `$ distributions for the ST and FT modes are shown in the upper and lower panels of FIG. 3, respectively. The WIMP signals are shown by the hatched histograms, the neutron backgrounds are shown in the filled histograms, and the total observable events are shown in the histograms with errorbars. We simply analyzed the asymmetry by comparing the number of events with $`\mathrm{cos}\gamma >\mathrm{cos}\gamma _0`$ ($`N_\mathrm{L}`$) and that with $`\mathrm{cos}\gamma <\mathrm{cos}\gamma _0`$ ($`N_\mathrm{S}`$), where $`\gamma _0`$ was chosen so that $`N_\mathrm{L}`$ equals $`N_\mathrm{S}`$ for a flat background. The asymmetry-detection confidence level($`ADCL`$) was then defined as

$$ADCL=\frac{|N_\mathrm{L}N_\mathrm{S}|}{\sqrt{N_\mathrm{L}+N_\mathrm{S}}}.$$
(3)
We obtained 5.8 $`\sigma `$ and 13 $`\sigma `$ in the ST mode and the FT mode, respectively. Although the ST mode showed a lower confidence level than the FT mode, it was shown that the ST mode is still a very strong method compared to the conventional method to detect the positive signal of WIMPs. Because the direction of the WIMP-wind has a diurnal modulation which goes out of phase with a period of a year, it is unlikely that any background that is isotropically distributed or localized near the detector, would mimic this asymmetry. Event the background that have or could have fixed day/night circadian rhythm does not mimic the WIMP signal because WIMP-wind is expected to be blowing in the same direction during the day now and during the night in six months.

We calculated $`ADCL`$ values for several parameter sets in order to evaluate its dependence on the experimental parameters. The experimental parameters given in the TABLE 1, except for the gas pressure and the energy threshold are used for the calculation. The result is shown in FIG. 4. It is seen that the ST mode and the FT mode have the largest $`ADCL`$ values at 30 Torr and 20 Torr, respectively. In FIG.5, 3 $`\sigma `$ detection sensitivities to the $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ as a function of $`M_{\mathrm{WIMP}}`$ are shown. The parameters given in TABLE 1, except for $`M_{\mathrm{WIMP}}`$, $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$, and the exposure, were used for the calculation. Here, the 3 $`\sigma `$ detection level was defined as the smallest cross section for which we observe $`ADCL`$ = 3 $`\sigma `$. Even the ST method can reach the best sensitivity of the current experiments with a $`0.3\mathrm{m}^3\mathrm{year}`$ of exposure at Kamioka Observatory. With the FT method, it it possible to explore the MSSM prediction region via SD interactions with a sufficient exposure ($`30\mathrm{m}^3\times 10\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}`$). When we see any implication of the WIMP-wind, a precise study can be performed with the same technology by analyzing the shape of the recoil angle distribution. A prototype WIMP detector with a detection volume of 30 $`\times `$ 30 $`\times `$ 30 $`\mathrm{cm}^3`$ is now being manufactured. Since the fundamental manufacturing technology is already established, a large-volume detector ($`1\mathrm{m}^3`$) for underground measurements will soon be available.

## 4 Conclusions

In this paper, we have shown that $`\mu `$-TPC filled with $`\mathrm{CF}_4`$ gas is a promising device for WIMP-wind detection via SD interactions. By the Full-Tracking method with sufficient exposure, it is expected that we can not only detect the WIMP-wind, but can also precisely study the nature of WIMPs.

## Acknowledgements

This work is supported by a Grant-in-Aid in Scientific Research of the Japan Ministry of Education, Culture, Science, Sports, Technology; Ground Research Announcement for Space Utilization promoted by Japan Space Forum; and Grant-in-Aid for the 21st Century COE ”Center for Diversity and Universality in Physics”


format: text

text:
# HST/NICMOS Imaging Polarimetry of PPNs

## 1. Introduction

We observed 4 PPNs of SOLE type (Ueta et al. 2000) by means of imaging polarimetry using HST/NICMOS (NIC1 with POL0S, POL120S, and POL240S filters) during cycle 11 (GO 9377). Our main objective was to study the dust distribution in these optically thin PPN shells around the bright central star. Imaging polarimetry (e.g., Gledhill et al. 2001) is a very effective technique for this purpose because (1) bright stellar emission (unpolarized component) is automatically removed when the linearly polarized component (i.e., dust-scattered star light of the reflection nebula) is computed from the Stokes $`Q`$ and $`U`$ parameters, (2) the dust distribution in the reflection nebula can be probed at high resolution as opposed to the mid-IR diffraction-limit resolution, and (3) dust properties (e.g., the size distribution) can also be investigated.

## 2. Preliminary Results

Figure 1 shows imaging polarimetry data for IRAS 07134+1005 (left) and IRAS 06530-0213 (right). The IRAS 07134+1005 data illustrate the distribution and strength of polarization. The four frames display, from top left to bottom right, the total intensity (the Stokes $`I`$), the polarized intensity ($`I_{\mathrm{pol}}=\sqrt{Q^2+U^2}`$), the degree of polarization ($`I_{\mathrm{pol}}/I`$), and the 10.3 $`\mu `$m Gemini image (Kwok et al. 2002) for comparison. The total and polarized intensity images look similar, suggesting that the nebula emission consists mostly of scattered star light as expected in this reflection nebula. The total intensity image confirms the equatorially enhanced (toroidal) structure of the shell as seen in the 10.3 $`\mu `$m image and also shows the structure within the reflection nebula in an unprecedented detail. Such detailed structure of the shell was not seen well in HST/WFPC2 data due to the invasively bright central star and its enormous PSF spikes (Ueta et al. 2001).

In the IRAS 06530-0213 data, the distribution and angles of polarization vectors are displayed. The four frames shows the the polarized intensity, the distribution of polarization vectors of less than 30% ($`I_{\mathrm{pol}}/I<0.3`$) strengths, of 30 to 40% strengths, and of more than 40% strengths, from top left to bottom right. The vectors are distributed in a clear centro-symmetric pattern with respect to the central star. This confirms the optically thin nature of this PPN shell. It is also apparent that the stronger polarization is found near the edge of the shell. This is because dust-scattered light near the nebula edge tends to have near 90 scattering angle, and hence, tends to yield the highest degree of polarization. Our data yielded rather high degrees of polarization in our target PPNs. This may be due to a population of small dust grains residing in these shells. We are currently finalizing our data analysis, and the results would shed light on the PPN shell structure as well as the properties of dust grains present in PPNs.

### Acknowledgments.

Support for GO-9377 was provided by NASA through a grant from the STScI, which is operated by the AURA, Inc., under NASA contract NAS 5-26555. Ueta was also supported by IUAP P5/36 financed by the OSTC of the Belgian Federal State.

## References

Gledhill, T. M., et al. 2001, MNRAS, 322, 321

Kwok, S., Volk, K., & Hrivnak, B. J. 2002, ApJ, 573, 720

Ueta, T., Meixner, M., & Bobrowsky, M. 2000, ApJ, 528, 861


text:
# Binary Neutron Star Inspiral Search in LIGO S1

## 1 Introduction

In 2002, the LIGO Scientific Collaboration and the LIGO laboratory organized the first science data run (S1). The Collaboration took data for 17 days, between 23 August and 9 September, using the three LIGO detectors (two detectors in the Hanford Observatory and one detector in the Livingston Observatory), and the GEO detector in Hannover, Germany. The detectors had not achieved their aimed sensitivities, and were at different completion levels with respect to their configurations, as described in . However, the noise level of at least some of the detectors was low enough to make them sensitive to inspiraling binary neutron stars in the Galaxy, and even the Magellanic clouds, making the data taking and analysis effort worthwhile and competitive with previous searches that produced upper limits on the rate of inspiraling binary neutron sources in the Galaxy ().

The results of the analysis of S1 LIGO data looking for gravitational waves from binary neutron stars were described in detail in . This article reports on some of the details of the data analysis done in , on the results obtained, and on some of the lessons learned that will lead to improved methods in the analysis of present and future science runs. At the time of writing, the Collaboration is analyzing the data from another 2-month science run, and is preparing to take data again with another 60+ days long run at the end of 2003.

## 2 Detectors’ sensitivity to binary neutron star systems in S1

A typical amplitude spectral density of the LIGO detectors’ noise during S1, interpreted as strain, is shown in Fig. 1. We can translate the noise spectral density into an optimal range, defined as the maximum distance at which a binary neutron star system, if optimally oriented and located, would be detected with a signal-to-noise ratio of 8. In average, the optimal range for the Livingston detector during S1 was 176 kpc, and for the Hanford 4km detector, 46 kpc. The optimal range for the 2km Hanford detector was slightly worse than for the 4km detector, and the data quality was not as good; the range for the GEO detector for the same sources was significantly worse. We then decided to use only data from the LIGO 4km detectors, which we call L1 (Livingston detector) and H1 (Hanford detector).

Seismic noise at the Livingston Observatory limited operations during most weekdays, so the overall L1 duty cycle during S1 was 42% (170 hours of data), while the duty cycle of H1 was 58% (235 hours). The time when the two detectors H1 and L1 were in operation amounted to 116 hours, representing a duty cycle of only 28%. We decided to analyze all of the data available from L1 and H1, including the times when a single detector was available. While we were not going to be able to confirm through coincidence the candidate events found when only one detector was operating, in this way we have more statistics to use for upper limit analysis and for testing our data analysis method. In the future, with longer science runs and improved duty cycles, we will be able to achieve good upper limits on astrophysical rates using only times when more than one detector is in operation.

The scientific goal of the data analysis exercise was to set an upper limit on the rate of binary neutron star systems in the Galaxy. We did not observe any coincident event with signal-to-noise ratio (SNR) larger than 6.5 in both detectors; thus, we had no candidates for detection.

In order to get an upper limit on the event rate, we measure the efficiency of the detectors (and the search method) to the population we are sensitive to (Milky Way Galaxy, Large and Small Magellanic Clouds), and infer the upper limit rate in that population from the efficiency as a function of signal-to-noise $`ϵ(\rho )`$ and the amount of time analyzed $`T`$. The result is an observational, frequentist, upper limit of 170 events/year per Milky Way Equivalent Galaxy (MWEG), with 90% confidence, on the coalescence rate of binary systems in which each component has a mass in the range 1-3 M.

A rough estimate of what could be expected from the analysis shows that our result was close to expectations. At SNR=16 (the SNR of our strongest surviving candidate), L1 was, in average, sensitive to 90% of the sources in the Galaxy, while H1 was sensitive to 40% of the same population. L1 was operating for 170 hours in S1, while H1 provided an additional 119 hours operating without L1 in coincidence. The total of available hours was then 289 hours. The maximum efficiency expected was then $`(0.9\times 170+0.4\times 119)/289=0.69`$. An upper limit on the event rate for the Galaxy is then expected at the level of $`2.3/ϵT`$ 100 events/yr. This is considerably smaller than the actual upper limit obtained of 170 events/yr. This was due to several reasons: not all the data was used because of availability, calibration problems or poor data quality; there was some loss of efficiency due to requirements on matching templates and on consistency under coincidence; and we used the upper bound in our rate estimate given the uncertainties in the results. Of course, even an optimally calculated event rate would be far from other astrophysical estimates which suggest rates of $`10^5`$/yr for our Galaxy.

## 3 Data analysis method and results

We developed a pipeline to analyze the data that could lead to detection during times when the two detectors were in operation, but also analyzed the data when a single detector was in operation. If this pipeline had resulted in a strong candidate appearing in both detectors, consistent with a coincident signal, we would have followed up with additional investigations to determine its origin, either astrophysical or instrumental. In S1, our analysis did not find any such coincident event. In the absence of detection, we take several more steps to produce, from the list of candidates, an observational upper limit on the rate of events in the Galaxy. We describe now the different elements leading to our resulting rate of 170 events/yr/MWEG.

The pipeline and details on the methods used for each step are described in detail in ; here we summarize the steps and point out possible improvements to apply in the analysis of future data.

### 3.1 Candidate events

The pipeline has several parameters that were tuned to obtain the best upper limit possible. Since biases could be introduced by this tuning, we performed the optimization on a selected data set, called the “playground”. This playground, about 10% of the data in coincidence, was not used in the final analysis. The playground times were chosen by hand, trying to choose locked segments that were representative of the different data quality over the run. However, a fair representation was difficult to obtain with such an ad-hoc procedure, and for the next science run (S2, Jan 14-Apr 14 2003), an automatic procedure was implemented.

In the first step in our pipeline, we analyze the interferometer data from each detector using matched filtering, with a template bank chosen to guarantee coverage of the worst detector  (and thus producing overcoverage of the better detector). The matched filtering makes critical use of the detector data calibration: we lost about 25 hours, or 9% of available data, due to missing calibration information. We also lost 39 hours, or 13%, due to the granularity of the matched filtering jobs, which required blocks of time 256 seconds long, and thus could not use short segments, or times at the end of segments. We expect to improve both of these numbers with better practice (for calibration), and better instrument stability (for longer operating segments).

We used a single template bank for both detectors, for all times in S1. In future runs, we will use different banks, dynamically adapted to the noise in the detector in the segment of continuous lock being analyzed. When the SNR of a particular candidate is larger than 6.5, we calculate a quality-of-fit parameter (similar to a $`\chi ^2`$, but adapted to the limitations placed by a finite template bank), and tune a cut on this parameter. This cut, although very powerful, was found not to veto candidates that upon further inspection were happening during noisy times. We are now developing methods that will analyze the time before the candidate template starts, either looking for excitations of other templates (typical during noisy times), or looking for lack of consistency with a stationary noise background.

The second step in the pipeline applied two kind of instrumental vetoes to the surviving triggers: an epoch cut, and an instrumental veto for H1. The epoch cut was adopted by the two working groups looking for signals of short duration (bursts and inspiral sources), and eliminated times when the noise in certain frequency bands was well outside typical noise levels during S1, trended over 6 minute periods. This cut removed 8% of the L1 data, and 18% of the H1 data. In spite of the non-negligible amount of data eliminated, we did not have a satisfactory explanation on the reasons why the noise was excessively high at these times. This was in part due to the lack of time to diagnose the noise sources: the instruments are constantly changing, except during the data taking run itself. This makes the diagnostics of particular data features very difficult, especially if attempted long after the instrument has changed character. We hope the understanding of the instrument will be better when we reach a stationary (and satisfactory) noise level, but probably not earlier. The second instrumental veto was due to coupling of frequency noise glitches in H1, as seen in another interferometer channel simultaneously with triggers in the gravitational wave channel. We had similar efficient veto channels for L1 (which also had a much higher event rate of accidental triggers), but we decided not to use them because when injecting hardware signals to simulate gravitational waves, they also excited the auxiliary channels in ways that could produce “vetoes”. We expect that with better understanding of the physical nature of the vetoes, we will be able to find efficient, safe vetoes that will eliminate the candidate events produced by instrumental artifacts.

The third step involves only the candidate triggers in L1 surviving the template matching and instrumental vetoes applied, and which are strong enough to show up in H1 with SNR$`>`$6.5. To these events, we apply a coincidence veto if the are no consistent triggers in H1 and L1. However, there were no such events in the S1 data: only triggers in L1 which appeared closer than 51 kpc would have possibly appeared as candidates in H1; but the strongest candidate event we had during coincident times had an apparent distance of 68 kpc. During S1, we were logistically limited in the number of triggers generated, thus imposing a lower limit of SNR looked for as 6.5. Thus, the lower limit in SNR for the noisier detector sets up a much stricter SNR criterion in the less noisy detector. To overcome this problem without overloading the database with low SNR events, we plan to implement a hierarchical search, looking for low SNR in the less sensitive detector only around times when a large SNR event is found in the more sensitive detector.

The final step in the pipeline generating the list of candidates involved maximizing all surviving triggers over time and over the template bank, to improve the timing resolution, and the confidence in the SNR of the candidate. This pipeline analyzed a total of 236 hours.

### 3.2 Detector efficiency

To measure the efficiency of our pipeline, we injected in software simulated signals, following the sample population for spatial and mass distributions from a Milky Way population produced by the simulations of Ref. , with the spatial distribution described in Ref. . Since L1 was sensitive to sources at a distance slightly larger than our Galaxy, we added sources from the Large and Small Magellanic Clouds, treating them as points at their known distances and sky positions. These systems contributed about 11% and 2% of the rate of the simulated signals, respectively. The total number of injected signals totaled more than 5000. The efficiency is also measured as a function of SNR: at each SNR, we measure the number of simulated sources in our model at the corresponding distance or closer. We measured an efficiency of 80% for SNR$`>`$8, and 53% for SNR=15.9, that of our loudest surviving candidate.

### 3.3 Upper Limit Result

In principle, we would like our pipeline to produce some non-zero rate of candidate events as a function of threshold SNR, and compare that with an estimated accidental rate at the same SNR. The background rate can be measured by time-shift analysis, which guarantees that any coincidence is accidental. However, our pipeline did not result in any event for which coincidence timing criteria could be applied, making this method for estimating the background impossible. For single detector data, no such reliable estimate of the background can be performed: we just have a list of candidates which we cannot confirm or veto through coincidence. Assuming the true sources have a Poisson distribution with rate $`R`$, then the probability $`P`$ of observing a signal with SNR=$`\rho `$ in a time of observation $`T`$, with the efficiency of the detector and method as a function of SNR given by $`ϵ(\rho )`$, is $`P=1e^{RTϵ(\rho )}`$. Given a list of candidate events with maximum SNR=$`\rho _{}`$, we interpret that with confidence level $`P`$, the rate is not larger than $`ln(1P)/ϵ(\rho _{})T`$, or $`2.3/ϵ(\rho _{})T`$ for 90% confidence.

Given our largest SNR=15.9 surviving in the L1-only pipeline, our measured efficiency $`ϵ(\rho =15.9)=0.53`$ to the simulated population, and our runtime $`T`$=236 hr, we obtain a rate $`R161`$ events/year. If we express the efficiency to sources in our Galaxy, comprising 88.5% of the population considered, we obtain an upper limit on the rate of 140 events/year/Milky Way Equivalent Galaxy. We have +14%/-10% uncertainties in the efficiency, mostly due to calibration uncertainties, and 5% uncertainties in the population estimates; we choose to state the obtained upper limit as the higher bound in our rate estimate, or 170 events/year/MWEG.

## 4 Conclusions

For the first time, we used data from coincident interferometric detectors to look for signals from binary inspiraling neutron star systems. We did not find any signal with SNR$`>`$6.5 while the two LIGO detectors were in operation. We used the data to set upper limits on the possible rates of inspirals in our Galaxy. These are very important successful milestones in the first steps towards the data analysis of interferometric data, and we have learned many important lessons from the exercise which we plan to use in the future analysis of science runs.

## Acknowledgments

We gratefully acknowledge the LIGO project and the LIGO Science Collaboration, who made the first LIGO Scientific Run possible. This work was supported by NSF grant PHY-0135389.

## References


text:
# THE DENSITY PROFILE OF MASSIVE GALAXY CLUSTERS FROM WEAK LENSING

## 1 Introduction

While the concordance flat $`\mathrm{\Lambda }`$CDM model, in which the matter density is dominated by cold dark matter (CDM), provides a good fit to observed large scale-properties of the universe, there remain some possible small-scale problems for this model.

Numerical simulations of structure formation in a CDM model predict that the dark matter (DM) halos of $`L_{}`$ galaxies such as the Milky Way should contain a number of subhalos that exceed the observed number of satellite dwarf galaxies by 1-2 orders of magnitude (e.g. Klypin et al. 1999; Moore et al. 1999a). Strongly suppressed star formation in the subhalos could be a possible solution to this problem. Observations of anomalous flux ratios of strongly gravitationally lensed multiple quasar images (Kochanek & Dalal 2003) and observations of the dynamics of optically dark high-velocity gas clouds in the local group (Robishaw, Simon & Blitz 2002) appear to be qualitatively consistent with this proposed solution.

In addition, the simulations predict that DM halos have cuspy inner density profiles $`\rho (r)r^\alpha `$, with $`\alpha `$ somewhere in the range between $`1.0`$ (Navarro, Frenk & White 1997; hereafter NFW) and $`1.5`$ (Moore et al. 1999b). This appears to contradict the observed dynamics of DM-dominated low surface brightness galaxies which favour softer cores with $`\alpha =0.2\pm 0.2`$ (de Blok, Bosma, & McGaugh 2003). On the scales of galaxy clusters, some studies indicate shallower density profiles than those predicted from CDM simulations (Sand et al. 2003), while others give $`\alpha `$ values that are consistent with CDM predictions (Bautz & Arabadjis 2003).

Attempts have been made to solve these small-scale problems of CDM by proposing DM models that modify its behavior on small scales. Some examples of these are models in which the DM is self-interacting (Spergel & Steinhardt 2000), self-annihilating (Kaplinghat, Knox & Turner 2000), fluid (Peebles 2000; Arbey, Lesgourgues & Salati 2003), warm (e.g., Sommer-Larsen & Dolgov 2001), repulsive (Goodman 2000), fuzzy (Hu, Barkana & Gruzinov 2000), decaying (Cen 2001), is both self-interacting and warm (Hannestad & Scherrer 2000), acts as mirror matter (Mohapatra, Nussinov & Teplitz 2002) or has its gravitational interaction with baryonic matter suppressed on small scales (Piazza & Marioni 2003). Of these, the self-interacting DM model of Spergel & Steinhardt is the one which has been explored in most detail. Here, we put limits on this model by using weak gravitational lensing to measure the average density profile of an ensemble of massive galaxy clusters. Details of this work are given by Dahle, Hannestad & Sommer-Larsen (2003).

## 2 Constraints on the DM halo profile

Our data set is a subset of the weak gravitational lensing measurements of 38 X-ray luminous clusters presented by Dahle et al. (2002). This subset consists of 6 clusters at $`z=0.3`$ for which weak gravitational shear has been measured out to a projected radius of 3 $`h_{65}^1`$ Mpc. We fit the average observed radial shear profile to a “generalized NFW profile” on the form

$$\rho (r)=\frac{\delta _c\rho _c}{(r/r_s)^\alpha (1+(r/r_s))^{3\alpha }}.$$
(1)
For the model above, the characteristic density $`\delta _c`$ is

$$\delta _c=\frac{200}{3}\left[_0^1x^2(cx)^\alpha (1+cx)^{\alpha 3}𝑑x\right]^1.$$
(2)
This model has a concentration parameter $`c`$ defined by $`c=r_{200}/r_s`$, where $`\overline{\rho }(r_{200})=200\rho _\mathrm{c}`$ and $`\rho _\mathrm{c}`$ is the critical density. The outer slope at $`r>>r_s`$, $`\rho r^3`$, is chosen to be the same as the outer slope of simulated CDM halos, while the inner slope $`\alpha `$ is a free parameter. The result of our fit is given in Figure 1 which shows joint confidence limits for $`\alpha `$ and $`c_{\mathrm{vir}}`$ (defined as $`c_{\mathrm{vir}}=r_{\mathrm{vir}}/r_s`$, where $`r_{\mathrm{vir}}`$ is the virial radius of the halo). The CDM simulations predict a $`(z+1)^1`$ redshift-dependence and significant intrinsic scatter in the values of $`c_{\mathrm{vir}}`$ (Jing 2000; Bullock et al. 2001). A prediction for the value and scatter of $`c_{\mathrm{vir}}`$ for massive ($`M_{\mathrm{vir}}=2\times 10^{15}M_{}`$) clusters in a $`\mathrm{\Lambda }`$CDM universe is indicated in Figure 1.

Our data constrain $`\alpha `$ to be in the range $`0.9<\alpha <1.6`$ (68% CL), and $`\alpha <0.5`$ is excluded at the 95% level. We also find that the data are consistent with an isothermal sphere with a finite core, $`\rho (r)(r^2+r_c^2)^1`$, where $`r_c`$ is the core radius. For the case of self-interacting dark matter, our constraints on the core radius implies a self-interaction cross section $`\sigma _{}0.1\mathrm{cm}^2\mathrm{g}^1`$ (c.f. Yoshida et al. 2000, Meneghetti et al. 2001). This is at least an order of magnitude smaller than the cross section required to explain the observed rotation curves of low surface brightness galaxies (Davé et al. 2001).

## 3 Comparison with other results

Different recent studies give a wide range of values for $`\alpha `$. Sand et al. (2003; see also Sand, these proceedings) use a combination of strong lensing data and spectroscopic measurements of stellar dynamics in the central galaxy of three clusters which contain both radial and tangential arcs to find an average value $`\alpha =0.52\pm 0.05`$, but they also find evidence for a significant scatter $`\mathrm{\Delta }\alpha 0.3`$. On the other hand, Bautz & Arabadjis (2003) find $`1<\alpha <2`$ and Lewis, Buote & Stocke (2003) find $`\alpha =1.19\pm 0.04`$, based on Chandra observations of the X-ray luminous intracluster medium in four clusters and in one cluster, respectively. In contrast to our weak lensing study (which only probe the DM density profile at radii where the baryonic component is not dynamically dominant), these strong lensing and X-ray studies are not directly comparable to simulations that only contain collisionless CDM. The above results indicate that future observational studies should simultaneously take into account both the baryonic component in stars and in the X-ray luminous intracluster medium as well as the DM. Similarly, all these components must be properly modeled in numerical simulations, if the simulations are to be directly compared to cluster observations on small ($`10`$ kpc) scales. In any case, all the recent studies indicate that the core sizes of massive clusters are too small to be consistent with any self-interacting dark matter having a cross section large enough to explain the rotation curves of dwarf galaxies.

Like previous weak lensing studies (e.g., Clowe & Schneider 2001, Hoekstra et al. 2002), we are not able to strongly distinguish between the outer slope of an isothermal sphere, $`\rho r^2`$, and the NFW slope $`\rho r^3`$. However, in a recent work, Kneib et al. (2003) use a combination of weak and strong gravitational lensing data based on HST imaging of the cluster CL0024+17 to find an outer slope $`>2.4`$. Their data is adequately fit by a NFW profile with $`c=22_5^{+9}`$, significantly higher than typical observed concentration parameters of rich clusters (e.g., Hoekstra et al. 2002; Katgert, Biviano & Mazure 2003), which are generally consistent with CDM predictions (see also Fig. 1). However, Chandra X-ray data (Ota et al. 2003), as well as dynamical studies based on galaxy spectroscopy (Czoske et al. 2002), indicate that this is not a fully relaxed, spherically symmetric system. Weak lensing measurements of a representative sample of dynamically relaxed clusters out to even larger radii than we probe in our study should eventually settle the issue of the value of the outer slope.

## Acknowledgments

I thank my collaborators Steen Hannestad and Jesper Sommer-Larsen, and acknowledge support from The Reseach Council of Norway through a post-doctoral research fellowship.

## References


text:
# Hot Gas in Planetary Nebulae

## 1. Introduction

The formation of planetary nebulae (PNe) is a complex issue. The current fast stellar wind from the central star of a PN must interact with the previous slow AGB wind; it has thus been suggested that PNe are formed by interacting stellar winds (Kwok, Purton, & Fitzgerald 1978). The frequently observed multi-polar and point-symmetric structures of young PNe suggest that collimated outflows are responsible for shaping PNe in early evolutionary stages (Sahai & Trauger 1998). The fast stellar winds and high-velocity collimated outflows, upon impacting the ambient medium, will be shocked to a temperature of $`T=(3/16)\mu v^2/k`$, where $`\mu `$ is the average mass per particle, $`v`$ is the shock velocity, and $`k`$ is the Boltzmann constant. Hot gas at X-ray-emitting temperatures, i.e., $`1\times 10^6`$ K, can be generated for shock velocities greater than $``$300 km s-1.

As reviewed by Guerrero et al. (in this volume), X-ray observations have revealed hot gas in a number of young PNe such as NGC 6543 and NGC 7027, and established stringent upper limits on the hot gas content in evolved PNe such as NGC 7293 (the Helix). The detection of diffuse X-ray emission from PNe provides a strong support for the interacting-stellar-winds scenario of PN formation. This formation mechanism is in fact similar to that for circumstellar bubbles blown by Wolf-Rayet (WR) stars. Massive WR stars are descendants of red supergiants or luminous blue variables which lose mass via copious slow winds; the fast WR wind sweeps up the previous slow wind and forms a circumstellar bubble. Comparisons between PNe and WR bubbles will help us understand both types of objects. In the first part of this paper, we will discuss the physical properties of the 106 K hot gas in PN interiors and make comparisons with shocked stellar winds in WR bubbles.

In the second part of this paper, we will discuss hot gas at cooler temperatures, i.e., a few $`\times 10^5`$ K, in order to: (1) determine whether some evolved PNe are filled with such cooler hot gas, and (2) study the interface layer between the 106 K interior gas and the 104 K nebular shell in young PNe. We will present far-UV observations of the 105 K gas in PNe and use the results to make predictions about the presence or absence of 106 K hot gas in PN interiors.

## 2. $`10^6`$ K Hot Gas

Chandra and XMM-Newton observations of diffuse X-ray emission from PNe allow us to examine the spatial distribution and temperatures of the hot gas. Among the six PNe with diffuse X-ray emission, Mz 3 and NGC 6543 are the only two that are adequately resolved by the instrumental point spread function, and both of them show a limb-brightened X-ray morphology, indicating that the dense X-ray-emitting gas is concentrated near the inner walls of their nebular shells. This is similar to what we have observed in the WR bubbles NGC 6888 and S 308 (Chu et al. 2003; Gruendl et al. 2004); see Fig. 1 for an XMM-Newton EPIC image of S 308 and a comparison with an optical image. The large angular size of S 308 makes it possible to see the detailed relationship between the hot interior gas and the cooler nebular shell.

X-ray spectra of the hot gas in PN interiors (Guerrero et al. in this volume) show a narrow range of temperatures, $`23\times 10^6`$ K. These temperatures are much lower than the $`10^7`$ K gas detected in superbubbles around young clusters, where colliding stellar winds may have produced the hot gas and the high temperatures (Townsley et al. 2003). On the other hand, the temperatures of the X-ray-emitting gas in PNe are similar to the $`12\times 10^6`$ K gas detected in the WR bubbles NGC 6888 and S 308 (see their spectra in Fig. 1).

Both PNe and WR bubbles have been modeled hydrodynamically (Zhekov & Perinotto 1996, 1998; García-Segura et al. 1996a,b), assuming that the fast stellar wind is shock-heated to high temperatures and forms a contact discontinuity with the photoionized 104 K nebular shell. It is further assumed that thermal conduction and mass evaporation take place at the hot/cold gas interface, which lowers the temperature and raises the density of the hot interior gas. In fact, the hot interior gas is dominated by the evaporated nebular material.

As the X-ray emission is proportional to the square of the density, the integrated X-ray spectrum of the hot gas will be dominated by the emission from the densest gas which is at the lowest temperatures. Therefore, the low temperatures derived from model fits to the X-ray spectra are expected. The X-ray luminosities of the diffuse emission from PNe and WR bubbles are, however, at least an order of magnitude lower than those expected from models (Chu et al. 2001; Wrigge et al. 1994), implying that the amount of hot gas is significantly less than that predicted by models with thermal conduction.

Several possibilities have been suggested to reconcile the discrepancies between observation and theory (Arnaud et al. 1996; Soker & Kastner 2003): (1) the fast stellar wind may have evolved with time, with a lower velocity or mass loss rate in the past; (2) the fast stellar wind is anisotropic, with a bipolar structure; and (3) the thermal conduction at the hot/cold gas interface may be inhibited by the presence of magnetic fields. Unfortunately, diffuse X-ray emission has been detected in only six PN interiors; two PNe are not adequately resolved spatially, and another two are detected with less than 100 counts. The current X-ray observations are clearly inadequate to constrain models and to discern the above possibilities. High-quality X-ray observations of PNe are badly needed.

## 3. 105 K Hot Gas

Hot gas at 105 K is traditionally studied with spectral lines of highly ionized species, such as C IV, N V, and O VI. If these ions are produced by thermal collisions, C IV, N V, and O VI will be the dominant ionization stages for C, N, and O at temperatures of $`1\times 10^5`$, $`2\times 10^5`$, and $`3\times 10^5`$ K, respectively. As the ionization potentials of C III, N IV, and O V are 47.9, 77.5, and 113.9 eV, photoionization of these ions becomes significant in PNe with central stars at stellar effective temperatures greater than $``$ 35,000, 75,000, and 125,000 K, respectively. Therefore, O VI provides the best diagnostic for 105 K gas, as it is the least easily confused by photoionization.

The O VI $`\lambda \lambda `$1032,1037 lines can be observed with the recently launched Far Ultraviolet Spectroscopic Explorer (FUSE). Four FUSE observing programs were carried out to search for 105 K gas in PN interiors. Some observations use spectra of the central stars to search for nebular O VI absorption, while others use spectra of the nebula (avoiding the central star) to search for O VI emission. These four programs are:

* O VI absorption observation of 10 PNe (PI: Gruendl; AO2)
* O VI emission observation of NGC 6543 (PI: Guerrero; AO3)
* O VI emission observation of NGC 7009 (PI: Iping; AO3)
* O VI emission observation of 9 PNe (PI: Guerrero; AO4)

These observations have produced many exciting results. First of all, FUSE observations detected O VI in NGC 6543 and NGC 7009, two PNe with known diffuse X-ray emission. The top panels of Fig. 2 display FUSE spectra of their central stars; both show pronounced P Cygni profiles in the stellar O VI lines, indicating the presence of fast stellar wind. NGC 6543 shows a stellar spectrum with nebular O VI absorption at a velocity blue-shifted from the nebular systemic velocity, indicating that it is associated with the interface layer on the approaching side of the nebular shell. NGC 7009 does not show obvious nebular O VI absorption in its stellar spectrum. The middle panels of Fig. 2 display FUSE spectra of apertures centered inside the nebulae. Both nebulae show obvious O VI emission lines. The bottom panels of Fig. 2 display FUSE spectra of apertures centered outside the nebulae, where only airglow lines are seen. FUSE observations of NGC 6543 and NGC 7009 demonstrate clearly that the O VI emission from interfaces between the 106 K interior and the 104 K nebular shell can be detected. Quantitative analysis of the O VI emission and absorption lines, complicated by the interstellar and circumstellar H2 absorption along the line of sight, is in progress.

FUSE observations of the Helix Nebula (Fig. 3), a PN confirmed to lack 106 K gas, show contrasting spectra in both the central star and the nebula. The spectrum of the central star shows stellar O VI absorption, but not P Cygni profiles, indicating a lack of fast stellar wind, which is consistent with the conclusion derived from IUE observations (Cerruti-Sola & Perinotto 1985). Neither O VI absorption nor O VI emission from the nebula is detected. It is thus reasonable to expect that a lack of both stellar P Cygni profile in O VI and nebular O VI is a good indication of an absence of hot gas in the PN interior. NGC 3132 is an awarded Chandra target in Cycle 5. The FUSE observations of its central star and the nebula, displayed in Fig. 3, show neither stellar P Cygni O VI nor nebular O VI. We predict that no 106 K hot gas exists in NGC 3132 and that Chandra observations of NGC 3132 will not detect any diffuse X-ray emission.

## 4. Future Work

X-ray observations have shown that hot gas exists in some PN interiors, and this interior hot gas may play an important role in the formation and evolution of PNe. To understand the interacting stellar winds in PNe, high S/N and high spatial resolution X-ray observations of more PNe are needed. FUSE observations of O VI emission offer an effective diagnostic for the presence of hot gas in PNe, and should be used to guide the target selection for X-ray observations. Future models of interacting stellar winds need to consider the plasma micro-physics and dynamical mixing of nebular material, in addition to a critical assessment of the degree of thermal conduction. When high-quality X-ray spectra are available, they need to be modeled with considerations of non-equilibrium ionization and temperature structure in the hot gas.

## References

Arnaud, K., Borkowski, K. J., & Harrington, J. P. 1996, ApJ, 462, L75

Cerruti-Sola, S., & Perinotto M. 1985, ApJ, 291, 237

Chu, Y.-H., Guerrero, M. A., Gruendl, R. A., et al. 2003, ApJ, in press

Chu, Y.-H., Guerrero, M. A., Gruendl, R. A., et al. 2001, ApJ, 553, L69

García-Segura, G., Langer, N., & Mac Low, M.-M. 1996a, A&A, 316, 133

García-Segura, G., Mac Low, M.-M., & Langer, N. 1996b, A&A, 305, 229

Gruendl, R. A., Guerrero, M. A., & Chu, Y.-H. 2004, in preparation

Kwok, S., Purton, C. R., & Fitzgerald, P. M. 1978, ApJ, 219, L125

Sahai, R. & Trauger, J. T. 1998, AJ, 116, 1357

Soker, N., & Kastner, J. H. 2003, ApJ, 583, 368

Townsley, L. K., Feigelson, E. D., Montmerle, T., et al. 2003, ApJ, 593, 874

Wrigge, M., Wendker, H. J., & Witsotzki, L. 1994, A&A, 286, 219

Zhekov, S. A., & Perinotto, M. P. 1996, A&A, 309, 648

Zhekov, S. A., & Perinotto, M. P. 1998, A&A, 334, 239


text:
# Accurate spline solutions of the Dirac equation with parity-nonconserving potential

## I Introduction

Consideration of the parity nonconservation (PNC) effects in atoms provides an important verification of the Standard Electroweak Model in the low energy sector. The most important information provides the neutral *Cs* atom, because the experimental and the theoretical accuracy is the best in this case. The analysis performed in  and  has indicated a deviation of the measured weak charge value $`Q_w`$ from that predicted by the Standard Model by 2.5 standard deviations $`\sigma `$. Later this value was diminishing and growing up again in a series of works -. It was understood that the radiative corrections to the PNC amplitude play an important role. This question was investigated in - where both electron self-energy (SE) and vacuum polarization (VP) corrections were evaluated by different methods. However, in these calculations $`\alpha Z`$ expansion for SE ($`\alpha `$ is the fine-structure constant, $`Z`$ is the charge of the nucleus), Uehling approximation for VP or other approximations were employed. To our mind the direct accurate calculations of SE and VP without the use of any approximations are desirable. For this purpose we propose in this paper the accurate spline solutions of the Dirac equation with the weak-interaction PNC potential.

Another possible area of application of PNC spline solutions could be the spectra of the Highly-Charged Ions (HCI). Though no experimental results are available here up to now, the analysis of the PNC amplitude in the two-electron ions can also provide an important test of the Standard Model -. The calculation of correlation effects in two-electron ions is much easier than in neutral atoms, but QED corrections are more important.

The rigorous direct way to the calculation of SE and VP corrections to the PNC amplitude (see ) is rather cumbersome. In order to avoid the difficulties and receive the accurate results we propose to use another way: accurate solutions of the Dirac equation with the weak-interaction PNC potential. Having got the solutions of the Dirac equation with this potential we can obtain self-energy corrections, using the standard methods  \-  based on the spline approximation ,. In this work we solve the first part of the problem: the obtaining of the complete basis set of accurate solutions of the Dirac equation with the weak PNC potential. The main difficulty is the extreme smallness of the weak interactions in the atomic scale. To overcome this difficulty we change the weak interaction potential by scaling it in such a way that it becomes significant but still small enough to apply the perturbation theory. Then we can obtain the necessary PNC matrix elements by simple rescaling. Actually our numerical procedure appear to be so accurate that the scaling will become necessary only in few cases, for the most singular operators and the enhancement of the weak interaction potential will not exceed 10 times.

## II Dirac equation with PNC weak interaction

We employ atomic units $`m_e=e=\mathrm{}=1`$ throughout the paper. Consider stationary Dirac equation with the weak-interaction PNC potential

$$\left[c\stackrel{}{\alpha }\stackrel{}{p}+V(r)+\beta c^2+V_w(r)\right]\psi (\stackrel{}{r})=E\psi (\stackrel{}{r})$$
(1)
where $`\stackrel{}{\alpha }`$,$`\beta `$ are Dirac matrices, $`\stackrel{}{p}`$, $`E`$ are the electron momentum and energy, $`c137.03599976`$ is the speed of light, $`V(r)`$ is the Coulomb potential of atomic nucleus (pointlike or extended) and $`V_w`$ is the weak-interaction PNC potential. This potential (see, for example, ) looks like

$$V_w=\gamma _5W(r)=\gamma _5\frac{G_F}{2\sqrt{2}}Q_w\rho _{nuc}(r)$$
(2)
$$\gamma _5=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)\begin{array}{c}Q_w=N+Z(\mathrm{\hspace{0.17em}1}\mathrm{\hspace{0.17em}4}\mathrm{sin}^2\theta _w),\hfill \\ \rho _{nuc}(r)=\rho _0\left[\mathrm{\hspace{0.17em}1}+\mathrm{exp}\left[\left(rc\right)/a\right]\right]^1\hfill \end{array}$$
(3)
where $`\theta _w`$ is the Weinberg angle, $`\mathrm{sin}^2\theta _w0.2230`$ , $`N`$ is the number of neutrons in the nucleus, $`Z`$ is the number of protons, $`G_F`$ is the Fermi constant and $`\rho _{nuc}(r)`$ is the Fermi distribution for nuclear charge density. The eigenfunctions of the Dirac equation with PNC potential have no definite parity, therefore we assume

$$\psi (\stackrel{}{r})=\left(\begin{array}{c}\phi (\stackrel{}{r})\\ \chi (\stackrel{}{r})\end{array}\right)=\frac{1}{r}\left(\begin{array}{c}g_{njl}^1(r)\mathrm{\Omega }_{jlM}(\stackrel{}{n})+ig_{njl}^2(r)\mathrm{\Omega }_{j\overline{l}M}(\stackrel{}{n})\\ if_{njl}^1(r)\mathrm{\Omega }_{j\overline{l}M}(\stackrel{}{n})+f_{njl}^2(r)\mathrm{\Omega }_{jlM}(\stackrel{}{n})\end{array}\right)$$
(4)
where $`\phi (\stackrel{}{r})`$ and $`\chi (\stackrel{}{r})`$ are the upper and lower components of Dirac bispinor, $`g_{njl}^1(r)`$, $`g_{njl}^2(r)`$, $`f_{njl}^1(r)`$ and $`f_{njl}^2(r)`$ are the radial functions, $`\mathrm{\Omega }_{jlM}`$ and $`\mathrm{\Omega }_{j\overline{l}M}`$ are spherical spinors with opposite parity. The equations for the components $`\phi `$, $`\chi `$ are

$$\left(EV(r)c^2\right)\phi (\stackrel{}{r})c(\stackrel{}{\sigma }\stackrel{}{p})\chi (\stackrel{}{r})W(r)\chi (\stackrel{}{r})=\mathrm{\hspace{0.17em}0}$$
(5)
$$\left(EV(r)+c^2\right)\chi (\stackrel{}{r})c(\stackrel{}{\sigma }\stackrel{}{p})\phi (\stackrel{}{r})W(r)\phi (\stackrel{}{r})=\mathrm{\hspace{0.17em}0}$$
(6)
Inserting (4) in (5),(6) and separating out the angular variables we get

$$\left[\begin{array}{cccc}V(r)& 0& c\left(\frac{d}{dr}\frac{\kappa _{jl}}{r}\right)& W(r)\\ 0& V(r)& W(r)& c\left(\frac{d}{dr}+\frac{\kappa _{jl}}{r}\right)\\ c\left(\frac{d}{dr}+\frac{\kappa _{jl}}{r}\right)& W(r)& V(r)2c^2& 0\\ W(r)& c\left(\frac{d}{dr}\frac{\kappa _{jl}}{r}\right)& 0& V(r)2c^2\end{array}\right]\left[\begin{array}{c}g^1\\ g^2\\ f^1\\ f^2\end{array}\right]=\epsilon \left[\begin{array}{c}g^1\\ g^2\\ f^1\\ f^2\end{array}\right]$$
(7)
where $`\epsilon =Ec^2`$ and

$$\kappa _{jl}=\{\begin{array}{cc}l\hfill & \text{if }j=l\frac{1}{2}\hfill \\ (l+\mathrm{\hspace{0.17em}1})\hfill & \text{if }j=l+\frac{1}{2}\hfill \end{array}$$
(8)
This system of equations for the radial functions was first derived in . However in  only one certain solution of Eq (7) was obtained. Our finite goal requires the knowledge of the complete system of the solutions of Eq (7). For this purpose we use the B-spline approach ,.

## III Parity-nonconserving B-spline Approach

Using the Galerkin method (see  and ) we express the system of equations (7) in terms of the least action principle. We use MIT bag-model boundary conditions , employed also in . In case of PNC functions these conditions look like

$$\{\begin{array}{cc}g^1(R)=f^1(R)\hfill & \\ g^2(R)=f^2(R)\hfill & \end{array}\{\begin{array}{cc}g^1(\mathrm{\hspace{0.17em}0})=\mathrm{\hspace{0.17em}0}\hfill & \\ g^2(\mathrm{\hspace{0.17em}0})=\mathrm{\hspace{0.17em}0}.\hfill & \end{array}$$
(9)
Here $`R`$ is the size of a spherical box, where an atom or an ion is enclosed. We choose action functional $`S`$ so that condition $`\delta S=0`$ leads to Eqs (7) and (9). Then we present the functions $`g^1,g^2,f^1\text{and}f^2`$ in the form of the linear combination of B-splines. The condition $`\delta S=0`$ reduces to the system of $`4n\times 4n`$ symmetric generalized eigenvalue equations. For constructing the B-spline system we use the grid $`r_i=\rho _i^4/Z,i=0,\mathrm{},N_{int}`$, that was described in . In our calculations we use the number of the grid intervals $`N_{int}=200`$, the order of B-splines $`k=9`$ and the parameter, which defines the density of knots near zero (see ), $`h_{set}=0.0167`$.

In this paper we investigate the hydrogenlike ion $`{}_{}{}^{238}U_{}^{91+}`$. The box size was taken to be $`R1.0`$ a.u. The parameters of the Fermi distribution are $`c=7.136`$ fm and $`a=2.3/4\mathrm{ln}3`$ . As a first test of the accuracy of our B-spline approach, we set $`V(r)=Z/r`$, $`N_w=0`$ and compare the energy values from B-spline approach with low-lying energy levels given by the Sommerfeld formula. The results are presented in the Table I in a.u. Our accuracy of calculation of the energy $`10^{14}`$ is the limit of accuracy for the calculation with Fortran double precision type of variables.

Another check of the spline accuracy is given in Table II, where the matrix elements of the operators $`r^m(m=2,1,0,1,2)`$ are compared on the $`2s_{1/2}`$ and $`2p_{1/2}`$ spline and exact (point-like nucleus) wave functions without weak interaction. The results in Tables I and II show that all the matrix elements $`<2s_{1/2}|r^m|2s_{1/2}>`$ and
$`<2p_{1/2}|r^m|2p_{1/2}>`$ apart from the case $`m=2`$ can be evaluated with the same relative accuracy $`10^{14}10^{15}`$ as the low-lying energies. Only in case of $`m=2`$ the accuracy is of order $`10^710^8`$ because of the singularity of $`r^2`$ operator.

## IV Scaling and the numerical tests

As the next step we increase the PNC weak interaction by scaling $`V_w=N_w\gamma _5W(r)`$, where $`N_w`$ is the scaling factor, and investigate the energy spectrum by changing $`N_w`$. Here we use the Fermi distribution for the nuclear charge density.

The scaling problem looks as follows: to find a scaling parameter $`N_w`$ for which 1) the contribution of weak interaction to the matrix elements of our interest is larger than numerical inaccuracy, 2) weak interaction is still small enough for using the perturbation theory.

The Table III presents the PNC weak interaction contribution to the energy of the states for the different values of $`N_w`$. In case $`N_w0`$ the states do not have a certain parity and we use the notations $`n\stackrel{~}{s},n\stackrel{~}{p}`$ formally classifying the different states by their origin. The analysis of the results of calculations shows that the dependence $`E=E(N_w)`$ is approximately parabolic. One can also obtain this result from the perturbation theory (PT):

$$|n\stackrel{~}{s}=|ns+N_w\underset{mp}{}\frac{mp|V_w(r)|ns}{E_{ns}E_{mp}}|mp+\frac{1}{2}\left(N_w\right)^2\underset{m}{}\underset{l}{}\frac{ls|V_w|mpmp|V_w|ns}{(E_{mp}E_{ls})(E_{ns}E_{mp})}|ls$$
(10)
$$n\stackrel{~}{s}|\stackrel{~}{H}|n\stackrel{~}{s}E_{ns}=\left(N_w\right)^2\underset{m}{}\frac{\left|mp|V_w|ns\right|^2}{E_{ns}E_{mp}}$$
(11)
Here the evident equalities $`ms|V_w|ns=0`$ and $`mp|H(r)|ns=0`$ were used. The right-hand side of Eq (11) is also presented in Table III for $`N_w=1`$. The comparison shows that we can obtain the PNC weak interaction contribution to the energy with accuracy $`0.03\%`$ if we use scaling parameter $`N_w=10^2÷10^3`$. Here by the numerical inaccuracy we undestand the B-spline method inaccuracy of calculation of nonperturbative energy value presented in the last column of Table III. This leads to the poor accuracy of the weak-interaction PNC contributions to the energy, obtained with the low values of scaling parameter $`N_w=1÷10`$. For $`N_w>10^4`$ higher orders of perturbation theory become to be not negligible.

The second-order PNC weak interaction correction to the energy is of no physical interest and was investigated here only for the demonstration of the possibilities of the scaling procedure. This procedure helps to evaluate the second-order PNC weak interaction correction even though this correction is beyond the level of accuracy of the spline approximation.

We should emphasize that the matrix elements of our interest are the matrix elements with linear dependence on weak interaction. SE and VP corrections to the PNC amplitude are the objects of this type.

For testing our system of PNC B-spline functions we calculated the set of the matrix elements $`2\stackrel{~}{s}|r^m|2\stackrel{~}{s}`$ with the first non-zero order of perturbation expansion, linear in $`V_w`$. The following tests were employed:

$$2\stackrel{~}{s}|r^m|2\stackrel{~}{p}=\underset{n}{}\frac{2\stackrel{~}{s}|r^{m/2}|n\stackrel{~}{s}n\stackrel{~}{s}|r^{m/2}|2\stackrel{~}{p}}{E_{2\stackrel{~}{p}}E_{n\stackrel{~}{s}}}+\underset{n}{}\frac{2\stackrel{~}{p}|r^{m/2}|n\stackrel{~}{p}n\stackrel{~}{s}|r^{m/2}|2\stackrel{~}{p}}{E_{2\stackrel{~}{p}}E_{n\stackrel{~}{p}}}$$
(12)
$$2\stackrel{~}{s}|r^m|2\stackrel{~}{p}=\underset{n}{}\frac{2s|r^m|nsns|V_w|2p}{E_{2p}E_{ns}}+\underset{n}{}\frac{2s|V_w|npnp|r^m|2p}{E_{2s}E_{np}}$$
(13)
With these tests we check the completeness of the PNC B-spline spectrum (Eq (12)) and the absence of the higher order PNC contributions after the introduction of the scaling parameter $`N_w`$ i.e. the possibility of rescaling (Eq (13)). The results are given in Table IV. The behaviour of the matrix elements with $`m=2,1`$ is similar. The most suitable scaling parameters for the calculation of these matrix elements are $`N_w=1`$ (the absence of the scaling) with relative accuracy $`110^8`$ and $`N_w=10`$ with relative accuracy $`510^8`$. The test of completeness of the spectrum for these operators gives the relative accuracy $`110^{12}`$. The most interesting matrix elements are the ones with $`m=1,2`$ because of their singularity. For $`N_w=1`$ the relative inaccuracy for the calculation of both matrix elements is about $`110^6`$, for $`N_w=10`$ it is $`110^7`$. The checking of completeness of the basis for this operators give the relative accuracy $`110^7`$ in the worst case. Summarizing these tests we conclude that the best scaling parameter for the calculation of matrix element with linear dependence on weak interaction is $`N_w=10`$. The accuracy of this calculation is about $`110^7`$.

Finally, we want to illustrate our method by calculation of VP corrections in the Uehling approximation to the matrix element of PNC weak interaction on $`|2s_{1/2}`$ and $`|2p_{1/2}`$ functions. Such electroweak radiative corrections for the HCI were first obtained in  by solving Dirac equation with the Uehling potential

$$V_u=\frac{2\alpha Z}{3\pi r}_1^{\mathrm{}}e^{2ry/\alpha }\left(1+\frac{1}{2y^2}\right)\frac{\sqrt{y^21}}{y^2}𝑑y.$$
(14)
The scheme of calculation of these corrections via PT looks as follows

$`2s^{}|V_w|2p^{}2s|V_w|2p={\displaystyle \underset{n2}{}}{\displaystyle \frac{2s|V_u|nsns|V_w|2p}{E_{2s}E_{ns}}}+{\displaystyle \underset{n2}{}}{\displaystyle \frac{2s|V_w|npnp|V_u|2p}{E_{2p}E_{np}}}+`$ (15)
$`+{\displaystyle \underset{n2}{}}{\displaystyle \underset{l2}{}}\left({\displaystyle \frac{2s|V_u|nsns|V_w|lplp|V_u|2p}{(E_{2s}E_{ns})(E_{2p}E_{lp})}}+{\displaystyle \frac{2s|V_u|lsls|V_u|nsns|V_w|2p}{(E_{2s}E_{ns})(E_{2s}E_{ls})}}+{\displaystyle \frac{2s|V_w|npnp|V_u|lplp|V_w|2p}{(E_{2p}E_{np})(E_{2p}E_{lp})}}\right)`$
$`{\displaystyle \underset{n2}{}}{\displaystyle \frac{2s|V_w|npnp|V_u|2p2p|V_u|2p}{(E_{2p}E_{np})^2}}{\displaystyle \underset{n2}{}}{\displaystyle \frac{2s|V_u|2s2s|V_u|nsns|V_w|2p}{(E_{2s}E_{ns})^2}}`$
$`{\displaystyle \underset{n2}{}}{\displaystyle \frac{2s|V_w|2p}{2}}\left[{\displaystyle \frac{|2s|V_u|ns|^2}{(E_{2s}E_{ns})^2}}+{\displaystyle \frac{|2p|V_u|np|^2}{(E_{2p}E_{np})^2}}\right]`$
where the notation $`|2s^{}`$ and $`|2p^{}`$ means that this functions were obtained as a solutions of the Dirac equation with the Uehling potential. The right-hand side of Eq (15) presents PT expansion. We include the second order of PT because of its significance. Our results for the matrix elements in the left-hand side of Eq (15) are: $`2s|V_w|2p=4.16398610^6`$ eV and $`2s^{}|V_w|2p^{}=4.21956610^6`$ eV. These values are slightly different from the values obtained in  because of our use of Fermi charge density distribution and more modern data for the nuclear radius. In total the left-hand side of Eq (15) is equal to $`5,558010^8`$ eV. The calculation in the first order of PT in the right-hand side of Eq (15) gives $`5.505810^8`$ eV. The evaluation of the second order terms yields the value $`5,227710^{10}`$ eV and the sum of the first and second orders terms is $`5.558110^8`$ eV. Thus the relative discrepancy between the left and right-hand sides of Eq (15) is about $`210^5`$. This is the limit of accuracy of the calculation of the difference between two matrix elements $`2s|V_w|2p`$ and $`2s^{}|V_w|2p^{}`$ that were calculated with accuracy $`10^7`$.

For our purposes it is interesting to solve the Dirac equation with the PNC weak interaction and to calculate then the matrix element of the Uehling potential in a direct way and via PT:

$$2\stackrel{~}{s}|V_u|2\stackrel{~}{p}=\underset{n}{}\frac{2s|V_u|nsns|V_w|2p}{E_{2p}E_{ns}}+\underset{n}{}\frac{2s|V_w|npnp|V_u|2p}{E_{2s}E_{np}}$$
(16)
The results of comparison of the left- and right-hand sides of Eq (16) are given in Table IV. The agreement is about $`110^7`$ for $`N_w=1÷10`$.

It may be interesting also to compare the evaluation of the Uehling-PNC matrix element in a way given by Eq (15) and by the application of Eq (16). For this comparison we rearrange the terms in the right-hand side of Eq (16), separating out the large contribution, containing small denominator $`\mathrm{\Delta }E^{(2)}=E_{2s}E_{2p}`$. We remind that in our calculation we do not include the radiative corrections in the Dirac energy values. Therefore $`\mathrm{\Delta }E^20`$ only due to the use of the non-Coulomb field of the nucleus (Fermi distribution). Using the smallness of the differences $`\mathrm{\Delta }E^{(n)}=E_{ns}E_{np}`$ we replace $`E_{ns}`$ by $`E_{np}`$ and vice versa in the right-hand side of Eq (16). Then we arrive at the equality

$$2\stackrel{~}{s}|V_u|2\stackrel{~}{p}\frac{\left[2p|V_u|2p2s|V_u|2s\right]2s|V_w|2p}{E_{2s}E_{2p}}\underset{n2}{}\frac{2s|V_u|nsns|V_w|2p}{E_{2s}E_{ns}}+\underset{n2}{}\frac{2s|V_w|npnp|V_u|2p}{E_{2p}E_{np}},$$
(17)
the right-hand side of which coincides with the first two terms in the right-hand side of Eq (15). Evaluation of the left-hand side of Eq (17) gives the value $`5.505710^8`$ eV in a good agreement with the first-order PT result in Eq (15).

The second order of PT can be obtained by calculation of the following expression

$`{\displaystyle \underset{n2}{}}{\displaystyle \frac{2\stackrel{~}{s}|V_u|n\stackrel{~}{s}n\stackrel{~}{s}|V_u|2\stackrel{~}{p}}{E_{2\stackrel{~}{s}}E_{n\stackrel{~}{s}}}}+{\displaystyle \underset{n2}{}}{\displaystyle \frac{2\stackrel{~}{s}|V_u|n\stackrel{~}{p}n\stackrel{~}{p}|V_u|2\stackrel{~}{p}}{E_{2\stackrel{~}{s}}E_{n\stackrel{~}{p}}}}`$ (18)
$`{\displaystyle \frac{2s|V_w|2p}{E_{2s}E_{2p}}}{\displaystyle \underset{n2}{}}\left({\displaystyle \frac{2p|V_u|npnp|V_u|2p}{E_{2p}E_{np}}}{\displaystyle \frac{2s|V_u|nsns|V_u|2s}{E_{2s}E_{ns}}}\right)`$
$`{\displaystyle \underset{n2}{}}{\displaystyle \frac{2s|V_w|2p}{2}}\left[{\displaystyle \frac{|2s|V_u|ns|^2}{(E_{2s}E_{ns})^2}}+{\displaystyle \frac{|2p|V_u|np|^2}{(E_{2p}E_{np})^2}}\right]`$
For this expression we obtain the value $`5,223910^{10}`$ eV. Then the sum of the first and second orders terms is $`5.557910^8`$ eV. Then, the relative accuracy of calculation of the sum of the first and second order terms is the same as in the first method: $`210^5`$. Here we also lose two orders in accuracy due to subtraction of the leading terms in Eq (17). In any case we can state emphasize that the accuracy of calculation of any PNC matrix elements is of order $`10^7`$. Of course, one can lose the accuracy of calculation because of subtraction of two large values. The question of the accuracy of calculation of some composite expression should be solved separately in each case. As we understand, the calculation of the radiative corrections to the emission PNC amplitude should not contain any subtractions.

Concluding, we can say that the method developed here presents a powerful tool for the calculation of any corrections to the PNC matrix elements in atoms and ions. The recommended scaling parameter is $`N_w=1`$ or $`N_w=10`$ and the relative accuracy for the calculation of PNC matrix elements achieved in all cases is not less than $`10^7`$. In particular, this approach can be applied to the calculation of the radiative corrections to the emission PNC amplitude with high accuracy without any approximations ($`\alpha Z`$ expansion, Uehling potential approximation and other). Including in our scheme the Dirac-Hartree-Fock potential we can calculate these radiative corrections also for neutral atoms, e.g. for $`Cs`$. It will help to solve finally the problem of the possible deviation of the measured weak charge value $`Q_w`$ from that predicted by the Standard Model.

###### Acknowledgements.

The authors thank I.Goidenko and O.Yu.Andreev for the help with the B-spline codes. They are also grateful to G.Soff for drawing their attention to the paper . The financial support by the RFBR grant 02-02-16758 and by Minobrazovanje grant E02-3.1-7 is acknowledged.