text:
# The Stellar Content of Obscured Galactic Giant H II Regions III.: W31

## 1 INTRODUCTION

Blum et al. (1999, hereafter Paper I) and Blum et al. (2000, hereafter Paper II) presented near infrared imaging and spectroscopic observations of two optically obscured Galactic giant H II regions, W43 and W42, respectively. These observations revealed the massive star clusters at the center of the H II regions which had been previously discovered and studied at longer wavelengths. Our earlier work builds off the original success by Hanson et al. (1997) in elucidating the massive stellar content in M17 through detailed spectral classification of O–type stars in the $`K`$band (Hanson et al., 1996). The combination of imaging and spectroscopy at near infrared wavelengths is proving to be a powerful tool in exploring the youngest stages of stellar evolution and the birth environments of massive stars. The chief advantage of the near infrared is that it is long enough in wavelength to lessen (by upto a factor of ten in magnitudes) the effect of interstellar extinction over optical wavelengths and still short enough in wavelength to potentially probe the stellar photospheric features of massive stars. Thus, the near infrared is ideal to explore stellar nurseries to identify and elucidate the nature of their massive stars.

In the present paper, we continue our investigation of the high end of the Galactic stellar mass function and search for massive young stellar objects (YSOs) with observations toward W31. W31 is a well known star forming region which has been observed at radio, far infrared (FIR), and mid infrared (MIR) wavelengths. The W31 complex consists of three H II regions, G10.2–0.3, G10.3–0.1, and G10.6-0.4 (Wilson, 1972). The first, G10.2–0.3 is a giant H II region (GHII) which we define (following the suggestion of Dr. Robert Kennicutt, private communication) as one which produces at least $`10^{50}`$ Lyman continuum ($`=`$ Lyc) photons per second. This is about ten times the luminosity of the Orion nebula and roughly the number emitted from the hottest single O3-type star. As these stars are not found in isolation, there is an implication that a “giant” H II region contains some minimum of multiple O-type stars. The embedded stellar cluster in G10.2–0.3 (which we refer to as W31) is the subject of this paper. Recent FIR and radio observations as well as a summary of earlier long wavelength work are presented by Ghosh et al. (1989).

W31 is located at $`l,b=`$ 10.2, $`0.3`$ in the first Galactic quadrant. Its distance is not well known ($`>`$ 4.1 $`\pm `$ 0.9 kpc) owing to the classic radio recombination line velocity degeneracy and a non-definitive molecular absorption velocity (Wilson, 1972). We have modified the Wilson (1972) lower limit by taking the sun–to–Galactic center distance as $`R_{}`$ $`=`$ 8 kpc (Reid, 1993). We discuss the distance to W31 in §4 in light of our new spectro–photometric observations. Smith et al. (1978) put the Lyman continuum output of W31 at 2$`\times `$1050 s-1 (for a distance of 4.1 kpc, corrected to $`R_{}`$ $`=`$ 8 kpc).

## 2 OBSERVATIONS AND DATA REDUCTION

$`J`$ ($`\lambda 1.28`$ $`\mu `$m, $`\mathrm{\Delta }\lambda 0.3`$ $`\mu `$m), $`H`$ ($`\lambda 1.63`$ $`\mu `$m, $`\mathrm{\Delta }\lambda 0.3`$ $`\mu `$m), and $`K`$ ($`\lambda 2.19`$ $`\mu `$m, $`\mathrm{\Delta }\lambda 0.4`$ $`\mu `$m) images of W31 were obtained on the nights of 04 May 1999 and 22 May 2000 with the f/14 tip–tilt system on the Cerro Tololo Interamerican Observatory (CTIO) 4m Blanco telescope using the facility infrared imager OSIRIS. 111OSIRIS is a collaborative project between the Ohio State University and CTIO. OSIRIS was developed through NSF grants AST 9016112 and AST 9218449. Spectroscopic data for bright stars in W31 were obtained on the nights of 19, 21, and 22 May 2000 using the f/14 tip–tilt system at the Blanco telescope with OSIRIS. OSIRIS is described in the instrument manual found on the CTIO web pages (www.ctio.noao.edu); see also DePoy et al. (1993). The tip-tilt system is described by Pérez & Elston (1998). The tip–tilt system uses three piezo–electric actuators to move the secondary mirror at high frequency in a computer controlled feed–back loop which corrects the natural image centroid motion. OSIRIS employs 0.16′′ pixels in the mode used for all data described in this paper.

All basic data reduction was accomplished using IRAF222IRAF is distributed by the National Optical Astronomy Observatories.. Each image/spectrum was flat–fielded using dome flats and then sky subtracted using a median combined image of five to six frames. For W31 itself, independent sky frames were obtained 500 arcseconds south of the cluster. Standard stars used the median combination of the data for sky.

### 2.1 Images

The OSIRIS 1999 May images were obtained under variable photometric conditions with some thin cirrus and in $``$ 0.5′′ to 0.6′′ FWHM seeing (with the tip–tilt correction). Total exposure times on source were 504 s, 281 s, and 181 s at $`J`$, $`H`$, and $`K`$, respectively. The individual $`J`$, $`H`$, and $`K`$ frames were shifted and combined (Figure 1$`a`$, 1$`b`$, 1$`c`$), and these combined frames have point sources with FWHM of $``$ 0.51′′, 0.54′′, and 0.56′′ at $`J`$, $`H`$, and $`K`$, respectively. DoPHOT (Schecter et al., 1993) photometry was performed on the combined images. The flux calibration was accomplished using standard star 9172 from Persson et al. (1998) which is on the Las Campanas Observatory photometric system (LCO) in combination with short $`J`$, $`H`$, and $`K`$ images on W31 on 22 May 2000 under photometric conditions. The LCO standards are essentially on the CIT/CTIO system (Elias et al., 1982), though color transformations exist between the two systems for redder stars. No transformation exists between OSIRIS and either CTIO/CIT or LCO; the impact of this lack of transformation is further discussed below in §3.4.

The standard observations were made immediately after the W31 data were obtained and within 0.1 airmass of the airmass for W31. No corrections were applied for these small differences in airmass. Aperture corrections using 20 pixel radius apertures were used to put the instrumental magnitudes on a flux scale using between five and eight relatively uncrowded stars on the W31 images.

Uncertainties for the final $`J`$, $`H`$, and $`K`$ magnitudes include the formal DoPHOT error added in quadrature to the published error of the photometric standards and to the uncertainty of the aperture corrections (between the OSIRIS 1999 and 2000 data and the OSIRIS 2000 data and the standard star). The sum in quadrature of the aperture correction and standard star uncertainties is $`\pm `$ 0.027, $`\pm `$ 0.018, $`\pm `$ 0.021 mag in $`J`$, $`H`$, and $`K`$, respectively. The DoPHOT errors ranged from approximately $`\pm `$ 0.01 mag to an arbitrary cut–off of 0.2 mag (stars with larger errors were excluded from further analysis).

The flat–field illumination for the 1999 May images was not uniform. A smooth gradient with full range of about 10$`\%`$ was present. Corrections to the final photometry were made based on observations of a standard star taken over a 49 position grid covering the array.

### 2.2 Spectra

The $`K`$band spectra of six of the brightest stars in the center of W31 were obtained with a 0.48′′ wide slit (oriented EW) in $`<`$0.5′′ to 7′′ FWHM seeing and divided by the spectrum of HR 6953 (B9V), HR 6962 (A2V), or HR 7183 (B8V) to remove telluric absorption features. Total on source integration times were between 10 and 21 minutes for each W31 object and about 100 seconds for each telluric standard. Br$`\gamma `$ absorption in the HR stars was removed by eye by drawing a line across it between two continuum points. One dimensional spectra were obtained by extracting and summing the flux in $`\pm `$ 2 pixel apertures (0.64′′ wide). The extractions include background subtraction from apertures centered $`<`$1.0′′ on either side of the object.

An $`H`$band spectrum was obtained of the brightest cluster object, W31 #1. The slit width, telluric correction, seeing, and extraction were similar to those for the $`K`$band spectra. HR 6962 (A2V) was used to correct for telluric absorption. The intrinsic Br absorption in HR 6962 was removed by making smooth fits to the line profiles, and then dividing the original spectrum by the fit line spectrum. Regions outside the line profiles were set to unity in the fit line spectrum. The total on source integration time was 14 minutes.

The wavelength calibration was accomplished by measuring the positions of bright OH- lines from the $`H`$band sky spectrum (Olivia & Origlia, 1992). Lines in the W31 spectra are identified by their relative differences between one and another. The measured dispersion is 0.0002762 $`\mu `$m pix-1 in $`H`$ and 0.0003683 $`\mu `$m pix-1 in $`K`$. The spectral resolution at 1.65 and 2.2 $`\mu `$m is $`\lambda /\mathrm{\Delta }\lambda `$ 3000.

## 3 RESULTS

The $`J`$, $`H`$, and $`K`$band images shown in Figure 1$`a`$1$`c`$ have been combined into a false color image (Figure 2) with colors red, green, and blue representing $`K`$, $`H`$, and $`J`$, respectively. Figure 2 reveals a rich cluster of newly formed stars at the heart of W31. Some of these (see below) are shown to be massive O stars from their $`K`$band spectra. Spectroscopic and photometric evidence is also given below which indicates that young stellar objects (YSO) are present; thus the W31 cluster must be very young and its stars still undergoing their birth pangs.

The bright object (source “A”) near the center of Figure 2 is a foreground star. Our unpublished $`K`$band spectrum and its position in the color–color plot (see below) suggest it is an M dwarf star. We have identified this object (and several others) on the Digital Sky Survey 333 Based on photographic data obtained using The UK Schmidt Telescope. The UK Schmidt Telescope was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council, until 1988 June, and thereafter by the Anglo-Australian Observatory. Original plate material is copyright (c) the Royal Observatory Edinburgh and the Anglo-Australian Observatory. The plates were processed into the present compressed digital form with their permission. The Digitized Sky Survey was produced at the Space Telescope Science Institute under US Government grant NAG W-2166. red image. The USNO A2 catalog gives the coordinates of this object as $`\alpha `$ (2000) = 18$`{}_{}{}^{\mathrm{h}}09_{}^{\mathrm{m}}26.32^\mathrm{m}`$, $`\delta (2000)=20`$$`19^{}13.7^{\prime \prime }`$.

### 3.1 Color Magnitude and Color–Color Diagrams

The $`HK`$ vs. $`K`$ color–magnitude diagram (CMD) and $`JH`$ vs. $`HK`$ color color plot are shown in Figures 3 and 4. A foreground sequence is seen in Figure 3 at $`HK`$ $``$ 0.5 mag. To the red and lying at $`HK`$ between about 1.0 and 2.0 is the W31 main sequence. Some stars for which spectra are presented below are labeled. The relative colors of the reddest stars can be used to distinguish two sequences in Figure 4. There is a well populated sequence beginning near $`JH`$ vs. $`HK`$ $`=`$ 0, 0 which is taken to be the reddened main sequence. A second sequence lying to the red in $`HK`$ is populated by sources with a clear indication of excess emission and much larger extinction. We argue below that these are young stellar objects (YSOs; e.g. Lada & Adams (1992); see also the discussion in Paper II), though the precise evolutionary state is difficult to pin down from near infrared observations alone. One of these objects, W31 #1, was observed spectroscopically and is discussed further in §3.2 and §4.

### 3.2 Spectra

Spectra for each of four hot stars in the W31 cluster are presented in Figure 5. We name these objects W31 #2, #3, #4, and #5 as they are the second to fourth brightest sources at $`K`$ in Figure 3, excluding foreground sources (defined as objects with $`HK`$ $`<`$ 0.6 mag). The $`H`$ plus $`K`$ spectrum of W31 #1 (the brightest cluster source) is shown in Figure 6. All the spectra have been ratioed by an A–type star as described in §2 to correct for telluric absorption; the spectra in Figure 5 have been further normalized by a low order fit to the continuum. In the case of W31 #1, the spectrum has been put on a flux scale by multiplication by an appropriate blackbody (10000 K). The $`H`$ and $`K`$band segments are normalized based on the photometry of Figure 3. The spectrum in Figure 6 is not corrected for extinction. The upturn in the continuum at the end of the $`H`$band is probably due to vignetting in the OSIRIS camera optics; the continuum slope in this wavelength region is not reliable though the line profiles should not be substaintially afected relative to the local continuum. An analogous downturn may be present at the very end of the $`K`$band.

#### 3.2.1 O Stars

The spectra of sources W31 #2–5 may be compared to the $`K`$band spectroscopic standards presented by Hanson et al. (1996). The features of greatest import for classification are (vacuum wavelengths) the C IV triplet at 2.0705, 2.0796, 2.0842 $`\mu `$m (emission), the N III complex at 2.116 $`\mu `$m (emission), and He II at 2.1891 $`\mu `$m (absorption). The 2.0842 $`\mu `$m line of C IV is typically weak and seen only in very high signal to noise spectra (Hanson et al., 1996). The presence of N III and He II leaves no doubt that these are O–type stars, and the C IV emission firmly places each in the kO5–O6 subclass. This class is predominantly made up of stars with optical MK spectral types of O5 to O6, but several O4 stars presented by Hanson et al. (1996) also show C IV emission. None of the spectra in Figure 5 show He I 2.12 $`\mu `$m absorption which is present in later type O7-9 stars.

Generally strong absorption in Br$`\gamma `$ is expected for kO5–6 dwarf and giant stars and weak absorption or emission for super–giants, so the absorption wings seen in Figure 5 suggests the stars are not highly evolved. However, there is also residual emission seen in W31 #2, 4, and 5. The emission appears superposed on a broader absorption profile and might otherwise be suggestive of a circumstellar origin (e.g. in a disk, though there is no evidence of excess emission for these stars in Figure 4) were it not for the following. There are two problems with Br$`\gamma `$ which make it difficult to interpret and thus render any observations speculative. The first is the fact that Br$`\gamma `$ has to be corrected in the telluric standard (§2), and this is only accurate to five or ten percent for the method used here. Since the correction is for absorption in the standard and is typically under corrected (the core of the line is straight forward to identify, but the wings and continuum more difficult), the tendency is to reduce the intrinsic absorption in the object. It is also possible that any residual absorption in the telluric standard would have a different profile and position than the object. The second problem is that the nebular emission, which is subtracted by using background apertures, can vary over small scales. Even though the background is an average from either side of the object and in close spatial proximity to it, there is no guarantee that it is uniform or linearly varying across the source. It is interesting to note that we find no residual emission from He I 2.06 $`\mu `$m. Though Br$`\gamma `$ is often seen in emission in O supergiants (Hanson et al., 1996), the character of the emission is not like that shown in Figure 5 with absorption superposed upon emission. In summary, we can not tell if the residual emission in the O–stars is real (i.e. circumstellar) or an artifact. Higher spectral and spatial resolution observations will be required to make further progress.

#### 3.2.2 W31 #1

The spectrum of W31 #1 is shown in Figure 6. The spectrum rises more steeply to the red relative to the unratioed spectra of the O stars discussed in §3.2.1 (only the continuum divided O star spectra are shown here). This suggests a larger line-of-sight extinction and/or excess emission at redward wavelengths than for the O stars. There are no discernible photospheric absorption features present in either the $`H`$ or $`K`$ bands, though the signal-to-noise is relatively high, $`>`$70, in the fainter $`H`$band. Besides the lines discussed above for the $`K`$band, lines of He I (1.700 $`\mu `$m) and He II (1.692 $`\mu `$m and 1.572 $`\mu `$m) are expected for luminous hot stars (Blum et al., 1997; Hanson et al., 1998). Fe II emission is identified at (1.6878 $`\mu `$m) and perhaps also at 1.7419 $`\mu `$m; see Figure 7. There is weak He I 2.06 $`\mu `$m emission seen in Figure 6; however He I 1.700 $`\mu `$m is absent from the $`H`$band spectrum of W31 #1. These facts suggest that the line producing region seen in W31 #1 is not very hot.

### 3.3 Mid and Far Infrared Sources

Our position for the stellar cluster places it near ($``$ 38′′NW) the strong far infrared peak S9 presented by Ghosh et al. (1989) and to the SE of the even more luminous far infrared source S7. The far infrared beam was several arcminutes in diameter and the positional uncertainty quoted is $`<`$ 1. The position given by Frogel et al. (1977) for the mid infrared source toward G10.2–0.3 (29′′ aperture diameter) is approximately 38′′ east of the near infrared cluster center; Frogel et al. (1977) quote coordinates good to about 20′′. S9 and the mid infrared source are thus possibly coincident with each other and also with the stellar cluster. However, we believe the positional offsets should be taken at face value. It is more likely that the stellar cluster has cleared away much of the gas and dust in its immediate vicinity and that the hottest stars are exciting the material surrounding this cavity. This situation is similar to that for W43 (Paper I). The dust distribution is apparently non–uniform and the exciting cluster gives rise to the strong far infrared sources S9 and S7 to the SE and NW of the cluster, respectively.

### 3.4 Extinction

The results of the preceding sections may be used to estimate the total line–of–sight extinction to individual sources in the W31 cluster. However, as discussed in § 2, the OSIRIS colors are not on a standard system. The largest effect is expected to be in the $`J`$ filter (Persson et al., 1998; Elston et al., 1997, Paper II) such that the observed $`JH`$ and $`JK`$ colors are more red than CIT/CTIO. There is no evidence for a significant color term in $`HK`$ (Persson et al., 1998).

An estimate of the extent to which color–terms are important can be made by comparing the W31 colors to the dashed reddening lines in Figure 4. The intrinsic colors for O stars and red giants are taken from Koornneef (1983) and Frogel et al. (1978), respectively. The former colors have been put on the CIT/CTIO system (Paper II) by making minor ($`<`$1$`\%`$) corrections; the latter colors are on the CIT/CTIO system. The reddening lines themselves are derived using the interstellar extinction law of Mathis (1990) which is itself based on average relations from multiple photometric systems. The range of values for the near infrared extinction laws considered by Mathis result in a difference of about 0.15 mag for the lines shown in Figure 4 at $`HK=`$ 2.0 mag. There is a trend for the OSIRIS colors to be too red in $`JH`$ evident in Figure 4 which is in agreement with expectations as noted above. If we were to apply the same transformation to the $`J`$band colors as for the CTIO facility IR imager CIRIM (Elston et al., 1997) which has similar filters and detector as OSIRIS, then the effect would be for the reddened main sequence of the cluster to more closely match the lower reddening line in Figure 4 and for the stars with excess to lie even further below the line. The correction would be about $``$0.2 mag in $`JH`$ for the red stars like W31 #1 in Figure 4. This suggests that the magnitude of the undetermined color transformation to put OSIRIS magnitudes on the CIT/CTIO system is consistent with a negligible correction for $`HK`$.

We will therefore assume that the $`HK`$ colors of OSIRIS are the same as for the CIT/CTIO system and compute estimates of the extinction based on $`HK`$ alone and the Mathis interstellar extinction law. Table 1 lists the photometry and derived $`A_K`$ for the four O-type stars presented in this paper. The mean extinction toward these four objects is $`A_K=`$ 1.71 mag $`\pm `$ 0.12 mag where the uncertainty is the standard deviation in the mean. We add to this value (in quadrature) 0.15 mag due to the expected range in the power–law exponent of the Mathis (1990) extinction law ($`\pm `$ 0.1 in the exponent).

## 4 DISCUSSION

W31 hosts a complex mixture of newly formed stars, gas, and dust as Figure 2 attests. The spectra presented in §3 indicate that some of these newly formed stars are quite massive and also that massive stars are still forming, or at least have not yet shed their natal material and revealed their photospheres. In the following sections we discuss these spectra and the stellar and young stellar content in more detail, establishing the basic properties of the W31 cluster.

### 4.1 Distance

As mentioned in §1, the distance (and hence luminosity) of W31 has not been well known. Wilson (1972) discusses the radio recombination line velocity for W31. The Galactic rotation model is degenerate for this case, but the absorption due to H2CO in intervening clouds is significantly higher than the H recombination line velocity indicating the absorption is not local and that the H II region is not at the near distance. Interestingly, Georgelin & Georgelin (1976) placed W31 at the far distance of 18.7 kpc; it is the most distant GHII region in their famous plot (their Figure 11) of Galactic structure, well beyond the Galactic center. On the other hand, the H2CO absorption does not extend to the maximum velocity on the line–of–sight, so only a minimum, intermediate distance, $`D>`$ 4.1 $`\pm `$ 0.7 kpc can be estimated from this method. Even this limit is quite uncertain as the H2CO may be affected by non–circular velocities.

We can now constrain the distance to W31 using the spectroscopic and photometric results of §3. We compute distances assuming the O stars shown in Figure 5 are ZAMS or dwarf luminosity class (i.e. hydrogen burning). The results of the following sections firmly demonstrate that W31 is a very young cluster; hence, we will not consider the O stars to be evolved giants for which ages would be between about 2.5 and 5 Myr (Schaller et al., 1992). The near infrared spectra of the O–stars are also not consistent with those of supergiants (§3.2.1). The spectral type in each case is assumed to be O5.5V, consistent with the infrared type assigned above, and only the luminosity is varied. For the ZAMS case, the $`M_K`$ is taken from Paper II. For the dwarf case, the distance is determined using the $`M_V`$ given by Vacca et al. (1996) and $`VK`$ from Koornneef (1983). In all cases the distance estimate uses the photometry and $`A_K`$ given in Table 1. The distance estimates are shown in Table 1 and result in mean distances of 3.1 $`\pm `$ 0.3 kpc and 3.7 $`\pm `$ 0.4 kpc for the ZAMS and dwarf cases, respectively. The uncertainty quoted for the mean distance is the standard deviation in the mean of the individual distances added in quadrature to the uncertainty in $`A_K`$ ($``$ 250–350 pc).

The Wilson (1972) distance is a lower limit. Comparing this limit with the spectro–photometric estimates above, the radio limit is consistent (within one $`\sigma `$) with the estimates using the dwarf or ZAMS luminosity. As mentioned above, we argue below that W31 is very young based on the presence of YSOs. Thus, we do not expect the O–type stars to be highly evolved and we can firmly constrain the distance to be within 3.1 and 3.7 kpc with the expectation that it may be at or near the ZAMS distance. We will adopt a mean distance of 3.4 kpc. Then, the region shown in Figure 1$`a`$ is about 1.7 pc across, and that of the central cluster 0.6 pc.

### 4.2 Luminosity and Age

We can estimate the age of the O stars by placing them in the Hertzprung–Russell diagram (HRD). The average O star luminosity was set in the previous section by adopting the intrinsic luminosities consistent with the observed spectral types (the average of the ZAMS and dwarf cases) in order to derive the distance to the cluster. To place the individual stars in the HRD, we use this average distance (3.4 kpc), the bolometric correction vs. spectral type relationship given by Vacca et al. (1996), the photometry of Table 1, and the $`VK`$ for O–type stars presented by Koornneef (1983). The result is shown in Figure 8. The uncertainty in $`M_{\mathrm{bol}}`$ is due to an assumed distance uncertainty of $`\pm `$ 0.3 kpc and the uncertainty in $`A_K`$ and $`K`$ as given in Table 1. The $`T_{\mathrm{eff}}`$ uncertainty is taken as the range in $`T_{\mathrm{eff}}`$ for the O5–O6 spectral types as tabulated by Vacca et al. (1996). The data points are compared to the model isochrones of Schaller et al. (1992). Figure 8 indicates that the W31 O–type stars are quite young, less than about one million years with a mean age of perhaps 500 thousand years.

The Lyman continuum (LyC) output of W31 cited by Smith et al. (1978) is 2$`\times `$1050 sec-1 (for $`R_{}=8`$ kpc). For a distance of 3.4 kpc to W31, this number becomes 1.4$`\times `$1050 sec-1 (14 O star equivalents – Vacca et al. (1996)). The LyC output for a single O5.5V star, similar to what we observe, is 2.7$`\times `$1049 sec-1. Thus, five such stars can account for the total ionizing radiation production estimated from radio measurements. For the distance of 3.4 kpc derived in the previous section, the average $`M_{\mathrm{bol}}`$ of the O stars is about $``$9.1 mag (Figure 8). This is about 0.2 mag fainter than the O5.5V luminosity, so the corresponding LyC value is also about 20$`\%`$ lower. If the O stars are truly on the ZAMS, then the distance, and hence average luminosity, is an additional 20$`\%`$ lower than the O5.5V value (i.e. 0.4 mag total). There is thus room for $``$ 5$``$9 O5.5 stars; i.e. we have found four hot stars which account for about 50 to 80 $`\%`$ of the ionization indicated by radio emission (which does not account for destruction by dust or loss by leakage). Ghosh et al. (1989) and Woodward et al. (1984) presented 5 GHz maps of the W31 region. They suggested that many of the sources in their maps could be embedded young stars. Some of these objects may contribute to the ionization in W31, though in light of the smaller distance derived above, none of their sources is particularly luminous (their most luminous source overlapping with our images would be consistent with a O9 or B0 star). We further discuss below the Ghosh et al. (1989) sources in conjunction with near infrared excess emission objects presented in this paper.

### 4.3 Star Formation

In the midst of the complex gas and dust morphology evident in Figure 2, we have identified four luminous O–type stars. Figure 3 indicates there are more which are not yet spectroscopically confirmed. And more importantly, Figure 4 indicates there are a number of red objects which might be still dramatically affected by their natal dust clouds. These are the objects which appear offset to the red from the normal sequence of stars in Figure 4; see the object labeled #1 in the Figure. A clear excess emission (at least in $`K`$) is indicated by the offset. This excess is attributed to accretion luminosity in a circumstellar disk for lower mass YSOs (Meyer et al., 1997; Lada & Adams, 1992) or perhaps circumstellar envelopes in the case of Herbig Ae/Be (HAEBE) stars (Hartmann et al., 1993; Pezzuto et al., 1997). Lada & Adams (1992) identified lower mass objects in this region of the color–color plot as possible protostars. The protostars are distinguished from more evolved PMS stars by their location farther along the reddening line indicating extinction by more obscuring circumstellar material. Accounting for a foreground extinction as derived for the O–type stars in Table 1, then the YSO candidates in Figure 4 would lie near the lower end of the region identified by Lada & Adams (1992) as being populated by protostars. However, Lada & Adams (1992) did not attempt to account for foreground extinction in their YSO sample.

In Papers I & II, we identified similar YSO candidates in W43 and W42. These objects are likely massive stars in the process of shedding their birth material. Like the YSOs identified in M17 by Hanson et al. (1997), the objects in W43 and W42 may be of somewhat lower mass (based on their position in the CMD) than the more luminous spectroscopically identified O stars (perhaps late O–type or early B–type; see Hanson et al. (1997)). In M17, the more massive stars showed no evidence for circumstellar disks, but several later type stars did. An obvious question is whether the most massive stars formed disks at all, but the case in M17 is not conclusive because the stars with disk signatures are spatially segregated from the most massive O–type stars.

Several possibilities exist. All the massive stars are coeval, but the most massive either formed with no disks or the active disk phase was rapid and now beyond detection. Or, perhaps the later types simply formed later in a sequential mode. It is significant that the object W31 #1 is both spatially located in the heart of the cluster and more luminous at $`K`$ than any of the O–stars, unlike the suggestion of a spatial segregation of the stars with disks and without in M17.

In any event, the near infrared excess emission alone does not indicate the physical source of emission (e.g. from a disk or envelope), but it is a clear indication that this is a YSO. W31 #1 is more than 0.5 mag brighter at $`K`$ than any of the O stars listed in Table 1 indicating it is a massive YSO.

The spectrum of W31 #1 shown in Figure 6 shows no sign of photospheric hot–star features. Therefore, it does not allow us to determine whether it is a ZAMS late-O or early-B–star which has yet to clear away its remnant natal material (like a compact H II region) or perhaps a younger object just transitioning from an earlier proto–stellar phase. The Brackett series emission suggests an internal ionizing source (recall that background apertures, §2.2, have been used in an attempt to remove the nebular emission from the circumstellar emission). The connection between dense, hot molecular cores which might be massive protostars in a collapse or accretion phase and ultra compact H II (UCHII) regions which represent enshrouded but hydrogen burning objects (albeit very young) is the subject of intense recent study; see Kurtz et al. (2000) for a review. These objects probably form an evolutionary sequence (Kurtz et al., 2000) and the former are often sought by observations toward the latter. Figure 7 shows a detailed view of the $`H`$band spectrum of W31 #1. There is both Brackett emission and Fe II emission present. Such permitted lines are indicative of dense flows or winds and require large amounts of UV photons and dense gas ($`>`$106 cm-3), for example, in the massive evolved star Eta Carina (Hamann et al., 1994). The Fe II emission is narrow, and may also be consistent with dense material in a circumstellar disk.

Thus the emission–line spectrum of W31 #1 is consistent with the picture that a massive object lies hidden within a large amount of circumstellar material. This object would be sufficiently massive to produce a dense outflow/inflow and to produce sufficient ionizing photons to ionize the flow/wind. Massive protostars may remain “hidden” (i.e. they are not detected as UCHII) during the collapse/accretion phase because the mass infall rate is sufficient to arrest the development of an H II region around the ionizing protostar (Walmsley, 1995). Since we see evidence for ionized circumstellar material around W31 #1, this suggests it is not still in such an early proto–stellar phase. In this case we might expect W31 #1 to be associated with one of the radio sources reported by Ghosh et al. (1989); i.e. it would be an UCHII region. In Figure 9 we plot an overlay of the $`K`$band image of Figure 1$`a`$ with the positions of the Ghosh et al. (1989) radio sources (5 GHz sources from their Table 2). The Ghosh et al. (1989) 5 GHz beam was 7.5′′$`\times `$4.2′′ FWHM. As can be seen in the figure, there is a radio source near the position of W31 #1. This is radio source 23 from Ghosh et al. (1989), and it is about 4′′ from the near infrared position of W31 #1. Thus, W31 #1 is a candidate UCHII region.

W31 #1 is not the only YSO candidate in the cluster; in fact, it is one of five bright objects in the excess region of the color–color plot (Figure 4). Each of these objects (W31 #1, 9, 15, 26, 30) is in the bright end of the $`K`$ vs. $`HK`$ CMD, so is potentially massive (see Table 2). They all appear to be forming in the central region of the cluster where massive stars may be preferentially expected; see, for example, Hillenbrand & Hartmann (1998). These objects are indicated in Figure 1$`a`$ with open circles. The position of W31 #26 is consistent with the 5 GHz radio source 21 of Ghosh et al. (1989). This is the brightest of the Ghosh et al. sources in our field (and one of the brightest in their entire list).

In Figure 9, we plot the positions of the remaining Ghosh et al. radio sources as well as their radio sources 21 and 23 discussed above. We searched our near infrared source list for candidate matching objects to each of the remaining Ghosh sources. We have identified the reddest object within the radio FWHM as a possible candidate and show the results in Table 3. We have identified a potential radio counterpart in all cases except for Ghosh radio source 15. Though there is no spectroscopic information for these sources, we suggest that they may also be UCHII regions based on their infrared colors and likely conincidence with the 5 GHz radio sources. All the candidate counterparts identified in Table 3 are sufficiently red that they have no $`J`$band detections (except for radio sources 21 and 23 as discussed above). Since the radio sources are relatively bright at $`K`$, the lack of $`J`$band detections suggests they also exhibit excess emission, though we can’t prove it. Radio source 15 has no obvious counterpart on the $`K`$band image. If the radio source represents an embedded object, then it is still too buried to be seen at shorter wavelengths and may represent the most recent, on–going, star birth in the W31 cluster. The last three candidates in Table 3 may be too faint to be associated with the radio sources, as the implied $`K`$ magnitudes are inconsistent with an interpretation as early B–type stars. If so, the true sources may still be buried and invisible at near infrared wavelengths, for example, as are the “western” UCHII regions in W49 – (Smith et al., 2000).

We have tabulated the Lyman continuum luminosity for each of the Ghosh et al. (1989) sources in our field in Table 3. These values have been corrected for the distance derived above, and so are somewhat less luminous than reported by Ghosh et al. (1989). We also tabulate spectral types based on these values using the derived properties of dwarf OB stars from Vacca et al. (1996). Most are about B0. Could W31 #1 harbor such an embedded source? If we assume the excess emission from W31 #1 is due to reprocessing of the stellar radiation by a circumstellar disk, we can estimate its intrinsic stellar luminosity. Using the tabulation of Hillenbrand et al. (1992) for reprocessed radiation from a flat circumstellar disk, we find that $`\mathrm{\Delta }K`$ $`=`$ 3.71 mag and $`\mathrm{\Delta }(HK)`$ $`=`$ $`0.54`$ mag. This is the maximum luminosity that such a reprocessing disk may produce and assumes that the disk is viewed along the line of sight (it is seen along an axis perpendicular to the plane of the disk) and extends to the stellar surface (there is no hole in the disk).

Adjusting the observed magnitudes for this putative disk, we find that the $`A_K`$ would be 2.9 mag for the Mathis (1990) extinction law. This results in an intrinsic $`K`$ magnitude which is 1.71 mag fainter than the O stars in W31, equivalent to a B0.5 star (Paper II). This luminosity or spectral type estimate is consistent with the Ghosh et al. (1989) estimate from the radio continuum observations; see Table 9. Such an object does not have a very hard ionizing continuum and so might also explain the lack of neutral helium emission in Figure 6. Low ratios ($`<10\%`$) of He I to Br$`\gamma `$ are expected for such stars (Blum & Damineli, 1999). Recall that one of the other YSO candidates, W31 #26, is associated with the brightest radio source in our field, radio source 21. An inclined circumstellar disk would explain why its apparent $`K`$band excess is not as prevalent as that for W31 #1 (see Figure 3 and Table 9).

While illustrative, our model for the production of the near infrared excess in W31 #1 is not unique. We have not considered accretion luminosity contributing to the $`K`$band emission or dust destruction of the ionizing radiation in the UCHII region, for example. We plan to make longer wavelength observations in the mid–infrared to further explore the nature of this object by expanding the wavelength baseline to better constrain the source of its circumstellar emission. In general, we envision a situation where the O–stars and most of the YSO candidates are revealed and are ionizing the nebula in W31. The majority of the Ghosh et al. sources, however, may be less evolved, optically thick, UCHII regions.

It is interesting to further speculate on the nature of the star birth history in W31. If the Ghosh et al. sources are indeed associated with the near infrared sources, then their radio luminosities (Table 3) indicate they are generally less massive than the spectroscopically identified O–stars. We then have a situation similar to what may have occured in M17 (Hanson et al., 1997) where the most massive stars are revealed first with the later type O and B stars taking longer to clear away their natal material (or being slower to form). In the present case, there may also be on–going star birth of somewhat less massive hot stars (the buried radio source 15 of Ghosh et al. (1989)). A further similarity to M17 is that the most massive O stars appear not to show obvious disk signatures, and the implication is that the active disk phase for the most massive O stars is shorter than for lower mass O and B stars or that these objects don’t form disks at all. In the present case we can not yet firmly place the YSO W31 #1 on a mass scale, but its near infrared and radio properties are consistent with an early B–type star surrounded by circumstellar material.

## 5 Summary

We have presented initial, yet detailed, observations of the stellar and young stellar content embedded in the giant H II region W31. Our $`J`$, $`H`$, and $`K`$ images reveal a dense stellar cluster in W31, including young stellar objects. $`K`$band spectra are presented for five objects. Four of these are identified as early O–type by comparison to standard stars. The fifth and brightest cluster member at $`K`$ has near infrared colors which are similar to known young stellar objects. Strong excess $`K`$band emission places it to the red of the normal stellar sequence in the $`JH`$ vs. $`HK`$ color–color plot. We find that the position of this object, W31 #1, is consistent with a 5 GHz thermal radio source and suggest it therefore may be an ultra compact H II region, though perhaps one with an early B–type exciting star. It has an $`H`$ and $`K`$band spectrum devoid of photospheric lines, indicating it is still heavily enshrouded by circumstellar material. The Brackett lines are in emission. Its $`H`$band spectrum exhibits permitted Fe II line emission which we suggest indicates a dense circumstellar environment, perhaps a strong inflow, outflow or circumstellar disk.

We identified four other YSO candidates based on photometry in the $`JH`$ vs. $`HK`$ color–color plot. These stars (which may also be massive as indicated by their $`K`$band magnitudes) and W31 #1 are found in the heart of the W31 cluster among the O–type stars. We also identified a thermal radio source with one of these YSO candidates. Further, for the remaining 5 GHz sources in our field, we identified potential $`K`$band counterparts with offsets less than the size of the FWHM of the radio observations ($`<`$5′′). We suggest that these may be ultra compact H II regions representing the most recent massive star birth in the cluster.

Spectral types of the newly identified O–stars and the photometry presented herein firmly constrain the distance to W31 which was previously uncertain from radio observations. Based on the $`K`$band spectra and presence of YSOs, we argued that the O–stars can not be highly evolved stars (i.e. giants or supergiants). We then computed distances assuming the O–star spectra were consistent with ZAMS or dwarf luminosity classes by taking the appropriate range of intrinsic luminosities for O–stars. We find the distance to W31 must be less than about 3.7 kpc which is consistent with the previous lower limit from radio studies. Shorter distances (about three kpc) would be inferred for the O–stars (but still consistent with the radio data within the uncertainties) if they are precisely on the ZAMS. The inferred age of the O–stars is about half a million years or less.

PSC appreciates continuing support from the National Science Foundation. AD thanks PRONEX and FAPESP for support. We appreciate the comments of an anonymous referee which have helped to improve our paper. We are grateful to S. Strom for useful conversations regarding young stellar objects.


format: text

text:
# 1 Introduction

## 1 Introduction

Lepton number (L) conservation is one of the most obscure sides of the standard model(SM) not supported by an underlying principle and following from an accidental interplay between gauge symmetry and field content. Any deviation from the SM structure may introduce L non-conservation($`L/`$). Over the years the possibility of lepton number non-conservation has been attracting a great deal of theoretical and experimental efforts since any positive experimental $`L/`$ signal would request physics beyond the SM. In addition it would also show, that neutrinos are Majorana particles .

Recent neutrino oscillation experiments practically established the presence of non-zero neutrino masses, a fact that itself points to physics beyond the SM. However neutrino oscillations are not sensitive to the nature of neutrino masses: they can either be Majorana or Dirac masses leading to the same observables.

The principal question if neutrinos are Majorana or Dirac particles can be answered only by studying a lepton number violating processes since $`\mathrm{\Delta }L=2`$ is a generic tag of Majorana neutrinos. Various lepton-number violating processes have been discussed in the literature in this respect (for recent review see ). They offer the possibility of probing different entries of the Majorana neutrino mass matrix $`M_{ij}^{(\nu )}`$. Among them there are a few $`L/`$ nuclear processes having prospects for experimental searches: neutrinoless double beta decay ($`0\nu \beta \beta `$), muon to positron $`(\mu ^{},e^+)`$ conversion and, probably, muon to antimuon $`(\mu ^{},\mu ^+)`$ conversion . They probe $`M_{ee}^{(\nu )}`$, $`M_{\mu e}^{(\nu )}`$ and $`M_{\mu \mu }^{(\nu )}`$ matrix elements respectively. Currently the most sensitive experiment intended to distinguish the Majorona nature of neutrinos are those searching for neutrinoless $`0\nu \beta \beta `$-decay  The nuclear theory side  of this process have been significantly improved in the last decade  allowing reliable extraction of fundamental particle physics parameters from experimental data. Muon to positron nuclear conversion $`(\mu ^{},e^+)`$ is another $`L/`$ nuclear process with good experimental prospects.

The important role of muon as a test particle in the search for new physics beyond the standard model has been recognized long time ago. When negative muons penetrate into matter they can be trapped to atomic orbits. Then the bound muon can disappear either decaying into an electron and two neutrinos or being captured by the nucleus, i.e., due to ordinary muon capture. These two processes conserving both total lepton number and lepton flavors have been well studied both theoretically and experimentally. However, there are two other not yet observed channels of muon capture: The muon-electron ($`\mu ^{},e^{}`$) and muon-positron ($`\mu ^{},e^+`$) conversions in nuclei :

$`(A,Z)+\mu _b^{}`$ $``$ $`e^{}+(A,Z)^{},`$ (1)
$`(A,Z)+\mu _b^{}`$ $``$ $`e^++(A,Z2)^{}.`$
Apparently, the ($`\mu ^{},e^+`$) and ($`\mu ^{},e^{}`$) conversion processes violate lepton number L and lepton flavor $`L_f`$ conservation respectively. Additional differences between the $`(\mu ^{},e^{})`$ and $`(\mu ^{},e^+)`$ lie on the nuclear physics side. The first process can proceed on one nucleon of the participating nucleus while the second process involves two nucleons as dictated by charge conservation . Note also that the $`(\mu ^{},e^{})`$ conversion amplitude is quadratic and $`(\mu ^{},e^+)`$ amplitude linear in the neutrino mass. Thus the second process looks more sensitive to the light neutrino masses. The present experimental limit on $`(\mu ^{},e^+)`$ conversion branching ratio in 48Ti is

$`R_{\mu e^+}(Ti)={\displaystyle \frac{\mathrm{\Gamma }(\mu ^{}+^{48}Tie^++^{40}Ca)}{\mathrm{\Gamma }(\mu ^{}+^{48}Ti\nu _\mu +^{48}Sc)}}<4.3\times 10^{12}.`$ (2)
In the present paper we study the light Majorana neutrino mechanism for the $`(\mu ^{},𝐞^+)`$ conversion. Despite the previous rough estimates  indicate a very small branching ratio for this mode of $`(\mu ^{},𝐞^+)`$ conversion, by far below the experimental bound (2), we undertake a detailed study of this mode for several reasons. First, the nuclear theory of $`(\mu ^{},𝐞^+)`$ conversion is not yet sufficiently elaborated, as in the case of $`0\nu \beta \beta `$-decay, and requires further development. Second, $`(\mu ^{},𝐞^+)`$ conversion may receive contribution from other mechanisms offered by various models beyond the SM such as the R-parity violating supersymmetric models, the leptoquark extensions of the SM etc. Some of these mechanisms may involve the light neutrino exchange and, therefore, from the view point of nuclear structure calculations they resemble the ordinary light neutrino mechanism. Thus our present study can be viewed as a first step towards a more general description of $`(\mu ^{},𝐞^+)`$ conversion including all the possible mechanisms.

Below we develop a detailed nuclear structure theory for the light neutrino exchange mechanism of this process on the basis of the nuclear proton-neutron renormalized Quasiparticle Random Phase Approximation (pn-QRPA) wave functions . We perform a realistic calculation of the width of this process for the nuclear target $`{}_{}{}^{48}Ti`$ using limits on neutrino masses and mixings from neutrino oscillation phenomenology. A comparison with the previous estimations of $`R_{\mu e^+}`$ will be also presented .

The paper is organized as follows. The possible values of Majorana neutrino masses and mixings are discussed in sect. 2. The amplitude and width of $`(\mu ^{},e^+)`$ conversion are derived in sect. 3. The details of the calculation for the case of $`(\mu ^{},e^+)`$ conversion in $`{}_{}{}^{48}Ti`$ and our results are given in sect. 4. In sect. 5 we summarize our conclusions.

## 2 Majorana neutrino mass matrix

The finite masses of neutrinos are tightly related to the problem of lepton flavor/number violation. The Dirac, Majorona and Dirac-Majorana neutrino mass terms in the Lagrangian offer different neutrino mixing schemes and allow various lepton number/flavor violating processes . The favored neutrino mixing schemes has to accommodate present neutrino phenomenology in a natural way, in particular, to answer the question of the smallness of neutrino masses compared to the charged lepton ones. The most prominent guiding principle in this problem is the see-saw mechanism which can be realized in various models beyond the SM. A generic neutrino mass term is given by the formula

$$^{D+M}=\underset{l,l^{}=e,\mu ,\tau }{}[\frac{1}{2}\overline{(\nu _{l^{}L}^{})^c}(M_L^M)_{l^{}l}\nu _{lL}^{}+\frac{1}{2}\overline{\nu _{l^{}R}^{}}(M_R^M)_{l^{}l}(\nu _{lR}^{})^c+\overline{\nu _{l^{}R}^{}}(M^D)_{l^{}l}\nu _{lL}^{}]+h.c.$$
(3)
The first two terms do not conserve the total lepton number $`L`$. Here, $`\nu _L^{}`$ and $`\nu _R^{}`$ are the weak doublet and singlet flavor eigenstates. The indices L and R refer to the left-handed and right-handed chirality states, respectively, and the superscript $`c`$ refers to the operation of charge conjugation. $`M_L^M`$ and $`M_R^M`$ are complex non-diagonal symmetrical 3x3 matrices. The flavor neutrino fields are superpositions of six Majorana fields $`\nu _i`$ with definite masses $`m_i`$. Yanagida, Gell-Mann, Ramond and Slansky suggested that the elements of $`M^D`$ and $`M_L^M`$ be comparable with the masses of charged leptons and the hypothetical scale of lepton number violation ($`M_{LNV}10^{12}GeV`$), respectively. Then by diagonalization of the Dirac-Majorana mass term one ends up with the three very light and three very heavy neutrino eigenstates. This is the celebrated see-saw mechanism. Enlarging the number of the right handed neutrino states $`\nu _R`$ one can introduce sterile light mass eigenstates which may play a certain role in the explanation of the neutrino oscillation data including the LSND results. However the active-sterile neutrino oscillations as a dominant channel seems to be disfavored according to a recent Super Kamiokande global analysis  and work in Tübingen .

Sticking to the three neutrino scenario one may try to reconstruct the corresponding mass matrix from the neutrino oscillation data. This requires certain assumptions on its structure or additional experimental data. Solar, atmospheric and LSND neutrino data give information on the neutrino mass square differences $`\mathrm{\Delta }m_{ij}^2`$ as well as on the mixing angles of the unitary matrix $`U`$  relating the weak $`\nu _{lL}^{}`$ and mass $`\nu _{iL}`$ neutrino eigenstates

$$\nu _{lL}^{}=\underset{i=1}{\overset{3}{}}U_{li}^{(\nu )}\nu _{iL}(l=e,\mu ,\tau ).$$
(4)
This information can be used to restore the neutrino mass matrix inverting its diagonalization as

$$^{ph}=U^{(\nu )}diag(m_1,m_2,m_3)U^{(\nu )T},$$
(5)
if additional assumptions about the overall mass scale as well as about the CP phases $`\zeta _{CP}^{(i)}`$ of the neutrino mass eigenstates are made. This matrix can be identified with the $`3\times 3`$ Majorana mass terms $`M_L^M`$ in (3) and used in various phenomenological applications, for instance, in analysis of lepton number violating processes. The recent literature contains many sophisticated studies made in this direction (see, for instance,  and references therein).

The elements of the Majorana mass matrix are related to the effective Majorana neutrino masses $`<m_\nu >_{\alpha \beta }`$ ($`\alpha ,\beta =e,\mu ,\tau `$) as

$$_{\alpha \beta }^{th}<m_\nu >_{\alpha \beta }=\underset{k}{\overset{3}{}}U_{\alpha k}^{(\nu )}U_{\beta k}^{(\nu )}\zeta _{CP}^{(k)}m_k,$$
(6)
neglecting mixing with heavy neutral states if they exist in the neutrino mass spectrum. Amplitudes of lepton number violating processes are proportional to the corresponding effective neutrino masses . Thus $`0\nu \beta \beta `$-decay amplitude is proportional to $`<m_\nu >_{ee}`$, the so called effective electron neutrino mass . From the currently most stringent lower limit on $`0\nu \beta \beta `$-decay half-life of $`{}_{}{}^{76}Ge`$ $`T_{1/2}^{0\nu }1.1\times 10^{25}years`$  one obtains $`<m_\nu >_{ee}<0.62eV`$ . The effective Majorana muon neutrino mass $`<m_\nu >_{\mu \mu }`$ is related with the light neutrino exchange modes of muonic analog of $`0\nu \beta \beta `$-decay , semipletonic decay of kaon $`K^+\pi ^{}\mu ^+\mu ^+`$  etc. $`<m_\nu >_{\mu e}`$ enters the amplitude of the $`(\mu ^{},e^+)`$ conversion  and of the kaon decay into a muon and a positron ($`K^+\pi ^{}\mu ^+e^+`$). Some other elements of $`^{th}`$ are associated with rare B-decays .

In Ref.  the maximal allowed values for the elements of $`^{ph}`$ have been deduced from solar, atmospheric, LSND data and restriction coming from $`0\nu \beta \beta `$-decay. The result is

$`\left(\begin{array}{ccc}<m_\nu >_{ee}& <m_\nu >_{e\mu }& <m_\nu >_{e\tau }\\ <m_\nu >_{\mu e}& <m_\nu >_{\mu \mu }& <m_\nu >_{\mu \tau }\\ <m_\nu >_{\tau e}& <m_\nu >_{\tau \mu }& <m_\nu >_{\tau \tau }\end{array}\right)`$ $``$ $`\left(\begin{array}{ccc}0.60& 0.97& 0.85\\ 0.97& 0.76& 0.80\\ 0.85& 0.80& 1.17\end{array}\right)\text{eV}.`$ (13)
Thus in our analysis of the light neutrino exchange mode of the $`(\mu ^{},e^+)`$ conversion process we shall assume $`|<m_\nu >_{\mu e}|0.97eV`$.

We shall also use a more conservative model independent estimate of the effective neutrino mass. Atmospheric and solar neutrino oscillation data show up $`\mathrm{\Delta }m^2<<(1eV)^2`$ suggesting that all the neutrino mass eigenstates are approximately degenerate at the 1 eV scale . This observation in combination with the tritium beta decay endpoint allows one to set upper bounds on masses of all the three neutrinos  $`m_{e,\mu ,\tau }3`$eV. Thus in the three neutrino scenario one derives $`|m_\nu _{ij}|9\text{eV}`$ for $`i,j=e,\mu ,\tau `$ .

## 3 The $`(\mu ^{},𝐞^+)`$ conversion mediated by light neutrinos

The process of $`(\mu ^{},𝐞^+)`$ conversion is very similar to the $`0\nu \beta \beta `$-decay. Both processes violate lepton number by two units and take place only if the neutrino is a Majorana particle with non-zero mass. However, there are other important differences: i) The available energies for these two processes differ considerably. In addition, the number of leptons in final states is different. These facts result in significantly different phase space integrals. ii) The emitted positron in $`(\mu ^{},𝐞^+)`$ conversion has large momentum and therefore the long-wave approximation is not valid. iii) As it will be shown below, in the case of light neutrino-exchange there is a singular behavior of the $`(\mu ^{},𝐞^+)`$ nuclear matrix element which is the additional source of the difficulties for the numerical integration. iv) In the case of the $`(\mu ^{},𝐞^+)`$ conversion there is a great number of nuclear final states. Nevertheless, the major contribution is here assumed to come from the transition to the ground state of the final nucleus.

We shall discuss the amplitude and width of $`(\mu ^{},e^+)`$ conversion in nuclei mediated by light Majorana neutrinos. This process is shown in Fig. 1. We concentrate only on the nuclear transition connecting the ground states of the initial and final nuclei, which is favored from the experimental point of view due to the minimal background. In this case the $`e^+`$ spectrum has a peak at the energy

$$E_{e^+}=m_\mu \epsilon _b(E_fE_i).$$
(14)
Here, $`m_\mu `$, $`\epsilon _b`$, $`E_i`$ and $`E_f`$ are the mass of muon, the muon atomic binding energy (for $`{}_{}{}^{48}Ti`$ $`\epsilon _b=1.45MeV`$), the energies of initial and final ground states, respectively. Latter on we that the kinetic energy of the final nucleus is negligible.

The weak interaction Hamiltonian in the neutrino mass eigenstate basis has the standard form

$$^{weak}(x)=U_{li}^{(\nu )}\frac{G_F}{\sqrt{2}}\left[\overline{l}_L(x)\gamma _\alpha (1+\gamma _5)\nu _{iL}(x)\right]j_\alpha (x)+h.c.(l=e,\mu ,\tau ),$$
(15)
where $`j_\alpha (x)`$ is the charged hadron current. The neutrino mixing matrix $`U_{li}^{(\nu )}`$ is defined in Eq. (4).

In second order of the weak-interaction we get for the $`(\mu ^{},e^+)`$ conversion the following matrix element

$`<f|S^{(2)}|i>`$ $`=`$ $`i\left({\displaystyle \frac{G_F}{\sqrt{2}}}\right)^2<m_\nu >_{\mu e}{\displaystyle \frac{1}{(2\pi )^{3/2}}}{\displaystyle \frac{1}{\sqrt{4E_\mu ^{}E_{e^+}}}}u^T(k_{e^+})C^1(1+\gamma _5)u(k_\mu ^{})\times `$ (16)
$`{\displaystyle \frac{2m_em_\mu }{4\pi m_\mu }}g_A^2_{<m_\nu >_{\mu e}}^\mathrm{\Phi }2\pi \delta (E_\mu ^{}+E_iE_fE_{e^+}),`$
where

$$_{<m_\nu >_{\mu e}}^\mathrm{\Phi }=\frac{M_F^\mathrm{\Phi }}{g_A^2}M_{GT}^\mathrm{\Phi },$$
(17)
with

$`M_F^\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{4\pi }{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d\stackrel{}{q}}{2q}}\times `$ (18)
$`{\displaystyle \underset{n}{}}({\displaystyle \frac{<0_i^+|_l\tau _l^+e^{i\stackrel{}{k}_{e^+}\stackrel{}{r}_l}e^{i\stackrel{}{q}\stackrel{}{r}_l}|n><n|_m\tau _m^+e^{i\stackrel{}{q}\stackrel{}{r}_m}\mathrm{\Phi }(r_m)|0_f^+>}{qE_\mu ^{}+E_nE_i+i\epsilon _n}}+`$
$`{\displaystyle \frac{<0_i^+|_m\tau _m^+e^{i\stackrel{}{q}\stackrel{}{r}_m}\mathrm{\Phi }(r_m)|n><n|_l\tau _l^+e^{i\stackrel{}{k}_{e^+}\stackrel{}{r}_l}e^{i\stackrel{}{q}\stackrel{}{r}_l}|0_f^+>}{q+E_{e^+}+E_nE_i+i\epsilon _n}}),`$
$`M_{GT}^\mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{4\pi }{(2\pi )^3}}{\displaystyle }{\displaystyle \frac{d\stackrel{}{q}}{2q}}\times `$ (19)
$`{\displaystyle \underset{n}{}}({\displaystyle \frac{<0_i^+|_l\tau _l^+\stackrel{}{\sigma }_le^{i\stackrel{}{k}_{e^+}\stackrel{}{r}_l}e^{i\stackrel{}{q}\stackrel{}{r}_l}|n><n|_m\tau _m^+\stackrel{}{\sigma }_me^{i\stackrel{}{q}\stackrel{}{r}_m}\mathrm{\Phi }(r_m)|0_f^+>}{qE_\mu ^{}+E_nE_i+i\epsilon _n}}+`$
$`{\displaystyle \frac{<0_i^+|_m\tau _m^+\stackrel{}{\sigma }_me^{i\stackrel{}{q}\stackrel{}{r}_m}\mathrm{\Phi }(r_m)|n><n|_l\tau _l^+\stackrel{}{\sigma }_le^{i\stackrel{}{k}_{e^+}\stackrel{}{r}_l}e^{i\stackrel{}{q}\stackrel{}{r}_l}|0_f^+>}{q+E_{e^+}+E_nE_i+i\epsilon _n}}).`$
Here, $`=r_0A^{1/3}`$ is the mean nuclear radius, with $`r_0=1.1fm`$ and $`m_e`$ is the mass of electron. $`\stackrel{}{r}_i`$ is a coordinate of the $`i`$th nucleon. $`E_\mu ^{}`$ ($`k_\mu ^{}`$) and $`E_{e^+}`$ ($`k_{e^+}`$) denote energies (four-momenta) of the bound muon ($`E_\mu ^{}=m_\mu \epsilon _b`$) and the emitted positron, respectively. $`E_n`$ and $`\epsilon _n`$ are respectively energy and width of the intermediate nuclear state. $`\mathrm{\Phi }(r)`$ is the radial part of bound muon in its orbit (see Appendix A).

In the derivation of nuclear matrix element $`_{<m_\nu >_{\mu e}}^\mathrm{\Phi }`$ we neglected the contribution from higher order terms of nucleon current (weak-magnetism, induced pseudoscalar coupling), which are expected to play a less important role. Following the analysis in Ref.  their consideration can reduce the value of $`_{<m_\nu >_{\mu e}}^\mathrm{\Phi }`$ by an amount of about $`20\%`$ by analogy to the $`0\nu \beta \beta `$-decay.

We have normalized the nuclear matrix element $`_{<m_\nu >_{\mu e}}^\mathrm{\Phi }`$ in the same way as usual for the corresponding $`0\nu \beta \beta `$-decay matrix element. We note that the denominators in the expressions for the $`(\mu ^{},e^+)`$ conversion and the $`0\nu \beta \beta `$-decay exhibit a different behavior. This is because the energy of the bound muon $`E_\mu ^{}`$ is large. The two denominators in Eq. (19) can be associated with the direct and the cross Feynman diagrams in Fig. 1. One notes that the value of $`(E_\mu ^{}+E_nE_i)`$ is negative. This fact implies that the widths of the nuclear states play an important role and that the imaginary part of the nuclear matrix element can be large. This point was not discussed in previous publications  and is one of the motivations of our $`(\mu ^{},e^+)`$ conversion calculation. We want to investigate if this singular behavior of the amplitude can lead to an enhancement of the $`(\mu ^{},e^+)`$ conversion branching ratio or not. In order to simplify the numerical calculations we complete the sum over virtual intermediate nuclear states by closure after replacing $`E_n`$, $`\epsilon _n`$ by some average values $`<E_n>`$, $`\epsilon `$, respectively:

$`{\displaystyle \underset{n}{}}{\displaystyle \frac{|n><n|}{qE_\mu ^{}+E_nE_i+i\epsilon _n}}`$ $`=`$ $`{\displaystyle \frac{1}{qE_\mu ^{}+<E_n>E_i+i\epsilon }},`$
$`{\displaystyle \underset{n}{}}{\displaystyle \frac{|n><n|}{q+E_{e^+}+E_nE_i+i\epsilon _n}}`$ $`=`$ $`{\displaystyle \frac{1}{q+E_{e^+}+<E_n>E_i+i\epsilon }}`$ (20)
Next we assume that the muon wave function varies very little inside the nuclear system, i.e., the following approximation is used

$$|_{<m_\nu >_{\mu e}}^\mathrm{\Phi }|^2=<\mathrm{\Phi }_\mu >^2|_{<m_\nu >_{\mu e}}|^2.$$
(21)
The explicit form of $`<\mathrm{\Phi }_\mu >^2`$ is given in Appendix B.

For the width of $`(\mu ^{},e^+)`$ conversion we obtain

$`\mathrm{\Gamma }_{<m_\nu >_{\mu e}}={\displaystyle \frac{1}{\pi }}E_{e^+}k_{e^+}F(Z2,E_{e^+})c_{\mu e}<\mathrm{\Phi }_\mu >^2|_{<m_\nu >_{\mu e}}|^2\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{m_e}}\right|^2,`$ (22)
where $`c_{\mu e}=2G_F^4[(m_em_\mu )/(4\pi m_\mu )]^2g_A^4`$ and $`k_{e^+}=|\stackrel{}{k}_{e^+}|`$. The nuclear matrix element $`_{<m_\nu >_{\mu e}}`$ can be decomposed into the contributions coming from direct and crossed Feynmann diagrams in Fig. 1 as

$$_{<m_\nu >_{\mu e}}=M^{dir.}+M^{cro.},$$
(23)
where

$`M^{dir.}`$ $`=`$ $`<0_i^+|{\displaystyle \underset{kl}{}}\tau _k^+\tau _l^+4\pi {\displaystyle \underset{\lambda }{}}(1)^\lambda \sqrt{2\lambda +1}j_\lambda (k_{e^+}R_{kl})\{Y_\lambda (\mathrm{\Omega }_r)Y_\lambda (\mathrm{\Omega }_R)\}_0\times `$
$`{\displaystyle \frac{}{\pi }}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{j_0(qr_{kl})j_\lambda (k_{e^+}r_{kl}/2)}{qE_\mu ^{}+E_nE_i+i\epsilon }}(\stackrel{}{\sigma _k}\stackrel{}{\sigma _l}f_A^2(q^2){\displaystyle \frac{f_V^2(q^2)}{g_A^2}})q𝑑q|0_f^+>`$
$`M^{cro.}`$ $`=`$ $`<0_i^+|{\displaystyle \underset{kl}{}}\tau _k^+\tau _l^+4\pi {\displaystyle \underset{\lambda }{}}(1)^\lambda \sqrt{2\lambda +1}j_\lambda (k_{e^+}R_{kl})\{Y_\lambda (\mathrm{\Omega }_r)Y_\lambda (\mathrm{\Omega }_R)\}_0\times `$ (24)
$`{\displaystyle \frac{}{\pi }}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{j_0(qr_{kl})j_\lambda (k_{e^+}r_{kl}/2)}{q+E_{e^+}+E_nE_i+i\epsilon }}(\stackrel{}{\sigma _k}\stackrel{}{\sigma _l}f_A^2(q^2){\displaystyle \frac{f_V^2(q^2)}{g_A^2}})q𝑑q|0_f^+>,`$
with

$`\stackrel{}{r}_{ij}=\stackrel{}{r}_i\stackrel{}{r}_j,r_{ij}=|\stackrel{}{r}_{ij}|,\stackrel{}{R}_{ij}=\stackrel{}{r}_i+\stackrel{}{r}_j,R_{ij}=|\stackrel{}{R}_{ij}|.`$ (25)
For the normalized nucleon form factors we use the conventional dipole form $`f_V(q^2)=1/(1+q^2/\mathrm{\Lambda }_V^2)^2`$ \[$`\mathrm{\Lambda }_V^2=0.71(GeV)^2`$\], $`f_A(q^2)=1/(1+q^2/\mathrm{\Lambda }_A^2)^2`$ \[$`\mathrm{\Lambda }_A=1.09GeV`$\].

## 4 Results and discussions

The nuclear matrix elements of the $`(\mu ^{},e^+)`$ conversion process have been calculated within the proton-neutron renormalized Quasiparticle Random Phase Approximation . The considered single-particle model space both for protons and neutrons have been as follows: The full $`03\mathrm{}\omega `$ shells plus $`2s_{1/2}`$, $`0g_{7/2}`$ and $`0g_{9/2}`$ levels. The single particle energies were obtained by using a Coulomb–corrected Woods–Saxon potential. Two-body G-matrix elements were calculated from the Bonn one-boson exchange potential within the Brueckner theory. The pairing interactions have been adjusted to fit the empirical pairing gaps . The particle-particle and particle-hole channels of the G-matrix interaction of the nuclear Hamiltonian $`H`$ are renormalized by introducing the parameters $`g_{pp}`$ and $`g_{ph}`$, respectively. The calculations have been performed for $`g_{ph}=1.0`$ and $`g_{ph}=0.8,1.0,1.2`$. The two-nucleon correlation effect has been considered in the same way as in Ref. .

In calculation of the $`(\mu ^{},e^+)`$ conversion nuclear matrix elements we have used the fact that the widths of the low lying nuclear states are negligible in comparison with their energies. Therefore we have carried out the calculation in the limit $`\epsilon 0`$ using the formula

$$\frac{1}{\alpha +i\epsilon }=𝒫\frac{1}{\alpha }i\pi \delta (\alpha ),$$
(26)
which allows one to separate the real and imaginary parts of the $`(\mu ^{},e^+)`$ conversion amplitude.

In Table 1 nuclear matrix elements of the light Majorana neutrino exchange mechanism of the $`(\mu ^{},e^+)`$ conversion in $`{}_{}{}^{48}Ti`$ are presented. The adopted value of $`<E_n>E_i`$ was 10 MeV. We have found that our results depend weakly on this average value of the nuclear states within the interval $`2MeV(<E_n>E_i)15MeV`$. However, they depend significantly on the details of nuclear model, in particular, on the renormalization of the particle-particle channel of the nuclear Hamiltonian, and on the two-nucleon short-range correlation effect (s.r.c.). A new feature of this $`(\mu ^{},e^+)`$ conversion calculation is that the imaginary part of $`_{<m_\nu >_{\mu e}}`$ is significant, i.e., can not be neglected. This fact was not noticed in the previous $`(\mu ^{},e^+)`$ calculations .

In the further analysis we shall consider the nuclear matrix element $`|_{<m_\nu >^{\mu e}}|`$ obtained for $`g_{pp}=1.0`$ by considering the two-nucleon short range correlations. It is interesting to compare its value with the value of $`0\nu \beta \beta `$-decay matrix elements for A=48 nuclear system. We have

$$|_{<m_\nu >_{\mu e}}|=2.45\times 10^2,|_{<m_\nu >_{ee}}|=0.82.$$
(27)
We see that the matrix element for the $`(\mu ^{},e^+)`$ conversion is strongly suppressed in comparison with $`0\nu \beta \beta `$-decay matrix element by about factor of 400. It is mostly due to the large momentum of the outgoing positron in the $`(\mu ^{},e^+)`$ conversion process.

One can compare also the width of $`(\mu ^{},e^+)`$ conversion in $`{}_{}{}^{48}Ti`$ with the width of $`0\nu \beta \beta `$-decay of $`{}_{}{}^{48}Ca`$. We get

$`{\displaystyle \frac{\mathrm{\Gamma }_{<m_\nu >_{\mu e^+}}}{\mathrm{\Gamma }_{<m_\nu >_{ee}}}}`$ $`=`$ $`{\displaystyle \frac{ln(2)}{G_{01}}}{\displaystyle \frac{1}{\pi }}E_{e^+}K_{e^+}F(Z2,E_{e^+})c_{\mu e}<\mathrm{\Phi }_\mu >^2{\displaystyle \frac{|_{<m_\nu >_{\mu e}}|^2}{|_{<m_\nu >_{ee}}|^2}}\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{<m_\nu >_{ee}}}\right|^2,`$ (28)
$`=`$ $`1.97\times 10^5{\displaystyle \frac{|_{<m_\nu >_{\mu e}}|^2}{|_{<m_\nu >_{ee}}|^2}}\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{<m_\nu >_{ee}}}\right|^2,`$
$`=`$ $`176.\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{<m_\nu >_{ee}}}\right|^2.`$
The width of the $`(\mu ^{},e^+)`$ conversion is enhanced mostly due to the larger available energy for this process. A comparison with the width of $`0\nu \beta \beta `$-decay show that it is disfavored by smaller coulombic factor F(Z,E) ($`0.623/1.8`$) and by significantly smaller value of associated nuclear matrix element ($`(0.0245/0.82)^2`$). If we assume the effective neutrino masses $`<m_\nu >_{\mu e}`$ and $`<m_\nu >_{ee}`$ to be comparable \[see Eq. (13)\] we find that $`\mathrm{\Gamma }_{<m_\nu >_{\mu e^+}}`$ is enhanced by a factor of about 200 in comparison with $`\mathrm{\Gamma }_{<m_\nu >_{ee}}`$. We have used $`G_{01}=8.031\times 10^{14}year^1`$ .

From the experimental point of view it is interesting to compare the $`(\mu ^{},e^+)`$ conversion width with the width of ordinary muon capture rate. We have

$`{\displaystyle \frac{\mathrm{\Gamma }_{<m_\nu >_{\mu e}}}{\mathrm{\Gamma }_\mu }}`$ $`=`$ $`2{\displaystyle \frac{E_{e^+}k_{e^+}}{m_\mu ^2}}{\displaystyle \frac{c_{\mu \mu }}{G_F^2}}{\displaystyle \frac{|_{<m_\nu >_{\mu e}}|^2F(Z2,E_{e^+})}{[G_V^2+3G_A^2+G_P^22G_AG_P]Zf(Z,A)}}\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{m_e}}\right|^2`$ (29)
$`=`$ $`2.24\times 10^{22}|_{<m_\nu >_{\mu e}}|^2\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{m_e}}\right|^2`$
$`=`$ $`1.34\times 10^{25}\left|{\displaystyle \frac{<m_\nu >_{\mu e}}{m_e}}\right|^2`$
If we use the prediction for $`<m_\nu >_{\mu e}`$ coming from neutrino oscillation phenomenology, i.e., $`<m_\nu >_{\mu e}0.97eV`$ and more conservative bound $`<m_\nu >_{\mu e}9eV`$ \[see Eq. (13)\], we end up with

$`{\displaystyle \frac{\mathrm{\Gamma }_{<m_\nu >_{\mu e}}}{\mathrm{\Gamma }_\mu }}=4.8\times 10^{37},4.5\times 10^{35},`$ (30)
This value is about ten orders of magnitude smaller as the estimated one for the $`(\mu ^{},e^+)`$ conversion in $`{}_{}{}^{32}S`$ by Doi et al. . There could be various reasons for this difference. First, the $`<m_\nu >_{\mu e}`$ nuclear matrix element calculated by Doi et al. contains contributions from all final nuclear states and not only from the $`0_{g.s.}^+0_{g.s.}^+`$ transition as in our case. It can be that the experimentally interesting $`g.s.g.s.`$ transition exhausts only a small part from all allowed nuclear transitions. Second, in the simplified calculation of Doi et al. nuclear matrix elements were evaluated by summing the final nuclear states with closure. It usually leads to overestimation of the results as we know from calculation of ordinary muon capture. Third, the nuclear matrix element of Ref.  have been evaluated by using the long–wave approximation. Our comparison of the $`(\mu ^{},e^+)`$ conversion and the $`0\nu \beta \beta `$-decay (long–wave approximation is used) matrix elements shows that it can lead to overestimation of $`_{<m_\nu >_{\mu e}}`$ by factor up to $`10^2`$. Fourth, the problem of the ground and short–range correlations have been not addressed in Ref. . We have found that $`_{<m_\nu >_{\mu e^+}}`$ matrix element for A=48 nuclear system is strongly suppressed by both of them. It is not clear whether this effect is due to the chosen target. However, we note that $`_{<m_\nu >_{\mu e^+}}`$ for the conversion in $`{}_{}{}^{48}Ti`$ consists of transition to the doubly closed shell nucleus $`{}_{}{}^{48}Ca`$, which, e.g., in the case of $`0\nu \beta \beta `$-decay are less favored. To clarify this issue, the calculations of the $`(\mu ^{},e^+)`$ conversion for other nuclear systems are necessary.

It is worthwhile to notice that our result is in relative good agreement with the calculations performed by Leontaris and Vergados  for the $`(\mu ^{},e^+)`$ conversion in $`{}_{}{}^{58}Ni`$. By using the same value for effective Majorana neutrino mass $`<m_\nu >_{\mu e}`$ as in this article the result of Ref.  corresponds to a branching ratio equal to $`3.2\times 10^{36}`$ relative to the total absorption of the muon for the ground state to ground state transition. The difference of about one order with our result for A=48 can be attributed to the nuclear physics aspect of the $`(\mu ^{},e^+)`$ conversion, i.e., to a given nuclear system and the chosen nuclear model. We also remark that in Ref.  the imaginary part of the $`(\mu ^{},e^+)`$ conversion amplitude was not considered.

## 5 Summary and outlook

In summary, the lepton number violating process of $`(\mu ^{},e^+)`$ conversion in nuclei have been studied. The light Majorana neutrino–exchange mechanism of this process have been considered. A detailed analysis of this mode of the $`(\mu ^{},e^+)`$ conversion have been performed. The first realistic calculation of the $`0_{g.s.}^+0_{g.s.}^+`$ channel of this process, which is most favored for experimentally studies, are presented. The relevant matrix elements for A=48 nuclear system have been calculated within the pn-RQRPA. The effects of the ground state and two-nucleon short-range correlations have been analyzed. It was found that by inclusion of them the value of $`(\mu ^{},e^+)`$ conversion matrix elements is strongly suppressed. We are the first, to our knowledge, showing that the imaginary part of the nuclear matrix element is large and should be taken into account. Further, a comparison of different relevant aspects with the $`0\nu \beta \beta `$-decay process are presented. It is shown that the width of $`(\mu ^{},e^+)`$ conversion is about by factor of 200 larger as that of $`0\nu \beta \beta `$-decay by assuming predictions for effective neutrino masses coming from neutrino oscillation phenomenology. Nevertheless, the studied neutrino exchange mode of lepton number violating $`(\mu ^{},e^+)`$ conversion is not suitable for experimental study being extremely small compared to the ordinary muon capture. This fact, however, does not disfavor further experimental study of the $`(\mu ^{},e^+)`$ conversion in nuclei as some other lepton number violating mechanisms, e.g., those coming from GUT’s and SUSY models, can dominate this process. Therefore, they should be carefully examine too.

Acknowledgments

We are grateful to I. Schmidt for useful comments and remarks. This work was supported in part by Fondecyt (Chile) under grant 8000017 and by RFBR (Russia) under grant 00-02-17587.

## 6 Appendix A

The bound muon wave function (1S-state) is

$$\mathrm{\Psi }_\mu (x)=\mathrm{\Phi }_\mu (\stackrel{}{x})e^{iE_\mu ^{}x_0}\frac{u_\mu ^s}{\sqrt{2E_\mu ^{}}},$$
(31)
where

$`\mathrm{\Phi }_\mu (\stackrel{}{x})`$ $`=`$ $`{\displaystyle \frac{Z^{3/2}}{(\pi a_\mu ^3)^{1/2}}}e^{Z|\stackrel{}{x}|/a_\mu },`$
$`u_\mu ^s`$ $`=`$ $`\left(\begin{array}{c}\chi ^s\\ 0\end{array}\right)\sqrt{2E_\mu ^{}}`$ (34)
with $`a_\mu =4\pi /(m_\mu e^2)`$ ($`a_\mu /a_em_e/m_\mu 5\times 10^3`$) $`m_\mu `$ is the reduced mass of muon nucleus system.

## 7 Appendix B

The width for the ordinary muon capture reaction $`\mu ^{}+(Z,A)\nu _\mu +(Z1,A)`$ can be written in the Primakoff form

$`\mathrm{\Gamma }_\mu ={\displaystyle \frac{1}{2\pi }}m_\mu ^2(G_F\mathrm{cos}\theta _c)^2<\mathrm{\Phi }_\mu >^2Z[G_V^2+3G_A^2+G_P^22G_AG_P]f(Z,A),`$ (35)
where muon average probability density over the nucleus is

$$<\mathrm{\Phi }_\mu >^2\frac{|\mathrm{\Phi }_\mu (\stackrel{}{x})|^2\rho (\stackrel{}{x})d^3x}{\rho (\stackrel{}{x})d^3x}.$$
(36)
$`\rho (\stackrel{}{x})`$ is the nuclear density. To a good approximation it has been found

$$<\mathrm{\Phi }_\mu >^2=\frac{\alpha ^3m_\mu ^3}{\pi }\frac{Z_{eff}^4}{Z},$$
(37)
i.e., the deviation from the behavior of the wave function at the origin has been taken into account by the effective proton number $`Z_{eff}`$. The values of this effective charge has been calculated for the nuclear systems of interest in Ref. . In particular, one finds $`Z_{eff}=17.6`$ for $`Z=22`$. The quadratic combination of the weak coupling constants is

$$[G_V^2+3G_A^2+G_P^22G_AG_P]5.9.$$
(38)
The function $`f(Z,A)`$ takes into account the two-nucleon correlations given by

$$f(A,Z)=10.03\frac{A}{2Z}+0.25(\frac{A}{2Z}1)+3.24(\frac{Z}{2A}\frac{1}{2}|\frac{1}{8Z}\frac{1}{4A}|).$$
(39)
This Pauli blocking factor for $`{}_{}{}^{48}Ti`$ takes the value $`f(22,48)=0.11`$.


text:
# Giant non-linearities accompanying electromagnetically induced transparency

## Abstract

We develop a fully quantum treatment of electromagnetically induced transparency (EIT) in a vapor of three-level $`\mathrm{\Lambda }`$-type atoms. Both the probe and coupling lasers with arbitrary intensities are quantized, and treated on the same footing. In addition to reproducing known results on ultraslow pulse propagation at the lowest order in the ratio of their Rabi frequencies, our treatment uncovers that the atomic medium with EIT exhibits giant Kerr as well as higher order non-linearities. Enhancement of many orders of magnitude is predicted for higher-order refractive-index coefficients.

PACS numbers: 42.50.Gy, 42.65.-k, 42.65.An

Introduction The discovery of electromagnetically induced transparency (EIT)  has led to the observation of new effects and development of new techniques in quantum optics. Recent examples include ultraslow light pulse propagation  and light storage  in atomic vapor, and atomic ground state cooling . It was proposed  that enhancement of efficiencies in nonlinear optical processes may be achieved by using EIT. A giant cross-Kerr non-linearity in EIT was suggested by Schmidt and Imamoǧlu , and has been indirectly measured in the experiment . We also note representative publications  utilizing EIT to study various optical non-linearities.

Essential to all these techniques is the phenomenon of EIT, in which coherent population trapping  is induced by atomic coherence and interference. As is well-known, a vapor of three-level $`\mathrm{\Lambda }`$-type atoms, initially in the ground state, when irradiated by two laser beams (the coupling and probe laser), exhibits EIT under certain conditions. Conventionally, both the coupling and probe lasers were treated as classical, external fields. In such treatments, the occurrence of EIT requires the coupling laser be much stronger than the probe laser. However, in the experiment  this condition was only marginally satisfied, which did not seem to have affected EIT. Moreover, with a quantum treatment of the probe laser, Fleischhauer and Lukin  have recently been able to predict the possibility of coherently controlling the propagation of light pulses via dark-state polaritons, which are formed through quantum entanglement of atomic and probe-photon states. This possibility is realized in the latest light storage experiments . This lesson teaches us that the quantum description of laser is more fundamental than the classical one, having advantages in uncovering new effects involving quantum nature of photons.

As a natural next step in this direction, in the present Letter we initiate a fully quantum treatment of EIT, namely we will deal with two quantized (coupling and probe) photon field modes interacting with three-level $`\mathrm{\Lambda }`$-type atoms. In this treatment the probe and coupling lasers are more or less on the same footing. This will have advantages in studying non-linearities in EIT, which requires the capability to deal with higher orders in the ratio of probe to coupling laser strengths (or the ratio of Rabi frequencies). Indeed, we will see that in addition to reproducing the well-known results on ultraslow light pulse propagation in atomic vapor at the lowest order, our treatment uncovers giant higher-order optical non-linearities in the atomic medium with EIT. They give rise to dramatic enhancement of Kerr as well as higher-order refractive-index coefficients.

The Model and Dressed States Let us consider a three-level atom, with energy levels $`E_1<E_3<E_2`$, interacting with two quantized laser fields, in the $`\mathrm{\Lambda }`$-type configuration (see Fig. 1): The lower two levels $`|1`$ and $`|3`$ are coupled to the upper level $`|2`$. Initially the atom is in the ground state $`|1`$. First applied is the coupling laser of frequency $`\omega _2=(E_2E_3)/\mathrm{}`$. It prepares necessary atomic coherence allowing the later applied probe laser of frequency $`\omega _1=(E_2E_1)/\mathrm{}\mathrm{\Delta }_1`$ to pass through the atomic medium without much absorption, even if the detuning $`\mathrm{\Delta }_1`$ is zero. This is the EIT  underlying the experimentally observed ultraslow pulse propagation in atomic vapor. In the following we will present a new treatment of EIT, in which both the probe and the coupling lasers are quantized.

With a unitary transformation $`U(t)=\mathrm{exp}(iH_0t/\mathrm{})`$, with $`H_0=_mE_m|mm|\mathrm{}[\mathrm{\Delta }_1(|22|+|33|)+\omega _1\widehat{a}_1^{}\widehat{a}_1+\omega _2\widehat{a}_2^{}\widehat{a}_2]`$, we adopt an interaction picture, in which the Hamiltonian of the atom-field system in the rotating-wave approximation remains time-independent:

$`\widehat{H}`$ $`=`$ $`\mathrm{}\mathrm{\Delta }_1|22|+\mathrm{}\mathrm{\Delta }_1|33|`$ (2)
$`+\mathrm{}(g_1\widehat{a}_1|21|+g_2\widehat{a}_2|23|+H.c.),`$
where $`|m`$ ($`m=1,2,3`$) are atomic states, $`g`$’s the coupling constants, and $`\widehat{a}_j`$ and $`\widehat{a}_j^{}`$ ($`j=1,2`$) the annihilation and creation operators of the probe and coupling laser modes. Assuming the detuning $`\mathrm{\Delta }_1`$ is small, we can develop a perturbation theory in which we retain only terms linear in $`\mathrm{\Delta }_1`$. The eigenstates of the Hamiltonian (2), also called the dressed states, have the following form:

$`|\varphi _{n_1,n_2}^{(i)}`$ $`=`$ $`a_i|1,n_1,n_2+b_i|2,n_11,n_2`$ (4)
$`+c_i|3,n_11,n_2+1,`$
where $`i=\pm ,0`$. The coefficients are given by

$`a_\pm `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_1}{\sqrt{2}\mathrm{\Omega }}}\left[1{\displaystyle \frac{\mathrm{\Omega }_1^2+4\mathrm{\Omega }_2^2}{2\mathrm{\Omega }^3}}\mathrm{\Delta }_1\right],`$ (5)
$`b_\pm `$ $`=`$ $`\pm {\displaystyle \frac{1}{\sqrt{2}}}\left[1\pm {\displaystyle \frac{\mathrm{\Omega }_1^2}{2\mathrm{\Omega }^3}}\mathrm{\Delta }_1\right],c_\pm ={\displaystyle \frac{\mathrm{\Omega }_2}{\sqrt{2}\mathrm{\Omega }}}\left[1\pm {\displaystyle \frac{3\mathrm{\Omega }_1^2}{2\mathrm{\Omega }^3}}\mathrm{\Delta }_1\right],`$ (6)
and

$`a_0={\displaystyle \frac{\mathrm{\Omega }_2}{\mathrm{\Omega }}},b_0={\displaystyle \frac{2\mathrm{\Omega }_1\mathrm{\Omega }_2}{\mathrm{\Omega }^3}}\mathrm{\Delta }_1,c_0={\displaystyle \frac{\mathrm{\Omega }_1}{\mathrm{\Omega }}},`$ (7)
with $`\mathrm{\Omega }_1=2g_1\sqrt{n_1}`$, $`\mathrm{\Omega }_2=2g_2\sqrt{n_2+1}`$ and $`\mathrm{\Omega }=\sqrt{\mathrm{\Omega }_1^2+\mathrm{\Omega }_2^2}`$. The energy eigenvalues of the dressed states (4) are

$`E_{n_1,n_2}^{(\pm )}={\displaystyle \frac{\mathrm{\Omega }_1^2+2\mathrm{\Omega }_2^2}{2\mathrm{\Omega }^2}}\mathrm{\Delta }_1\pm {\displaystyle \frac{\mathrm{\Omega }}{2}},E_{n_1,n_2}^{(0)}={\displaystyle \frac{\mathrm{\Omega }_1^2}{\mathrm{\Omega }^2}}\mathrm{\Delta }_1.`$ (8)
The complete set of dressed states for the system under consideration comprises the states $`|\varphi _{n_1,n_2}^{(i)}(i=\pm ,0)`$ for each $`n_1>0`$, and $`n_20`$ together with other two states $`|1,0,n_2`$ and $`|3,n_1,0`$ with zero eigenvalues. Knowing this set allows us to determine the time evolution of the atom-field system for any initial configuration.

Steady State and Susceptibility Assume that the atom is initially in the ground state, while the coupling and probe lasers in a coherent state $`|\alpha ,\beta `$, with $`\alpha `$ and $`\beta `$ supposed to be real for simplicity. Namely the initial state of the atom-field system is assumed to be

$$|\mathrm{\Psi }(0)=|1|\alpha ,\beta .$$
(9)
In the usual treatment with both lasers being classical external fields, one is concerned with a steady state of the atomic system. The counterpart of the steady state in our fully quantum treatment, according to the commonly used adiabatic hypothesis in quantum scattering theory and quantum field theory, is the state that evolves from the initial state (9) with the couplings $`g_1`$ and $`g_2`$ adiabatically turned on. Physically this is equivalent to having localized laser pulses before they enter the atomic vapor, with the pulse shape sufficiently smooth. In conformity to the Adiabatic Theorem, we need to identify a linear combination of the dressed states (4) that tends to the initial state (9) if we take the limits $`g_1,g_20`$ (or $`\mathrm{\Omega }_1,\mathrm{\Omega }_20`$). According to Eqs. (5-8), the ordering of the limits $`\mathrm{\Omega }_10`$ and $`\mathrm{\Omega }_20`$ is important. Corresponding to the actual conditions in which EIT is observed, the correct ordering is first $`\mathrm{\Omega }_10`$ and then $`\mathrm{\Omega }_20`$. In our interaction picture, this procedure selects the $`i=0`$ state in Eq. (4). Transforming it back to the Schrödinger picture, we identify the following state as the state that evolves adiabatically from the initial state:

$`|\mathrm{\Psi }(t)={\displaystyle \underset{n_1,n_2=0}{\overset{\mathrm{}}{}}}e^{(\alpha ^2+\beta ^2)/2}e^{iEt}{\displaystyle \frac{\alpha ^{n_1}\beta ^{n_2}}{\sqrt{n_1!n_2!}}}|\varphi _{n_1,n_2}^{(0)},`$ (10)
with $`E=E_1+\omega _1n_1+\omega _2n_2+E_{n_1,n_2}^{(0)}`$. This state can be viewed as the counterpart of the usual ”steady” state in our treatment. We see that the population of the upper level $`|2`$ is zero up to first order in detuning $`\mathrm{\Delta }_1`$. This means that there is no absorption, implying the phenomenon of EIT. We note that for EIT to occur the usual treatment requires the coupling laser be much stronger than the probe one , i.e. $`\mathrm{\Omega }_2\mathrm{\Omega }_1`$; then most atoms are populated in the ground state. Our formalism does not require this for EIT, so broader conditions are allowed: The coupling and probe lasers can be equally strong or even the probe laser is stronger. Then atoms can be in a more general superposition of two atomic states, which are entangled with the coupling and probe fields.

From Eq. (10) one can form the total density operator of the atom-field system. After taking trace over the field states, one gets the atomic reduced density operator. Then the Fourier component of the optical coherence of the ”steady” state at the frequency of the probe laser is found to be

$$\rho _{21}^A(\omega _1)=a_0(\overline{n}_\alpha ,\overline{n}_\beta )b_0(\overline{n}_\alpha ,\overline{n}_\beta ),$$
(11)
where $`\overline{n}_\alpha =\alpha ^2`$, and $`\overline{n}_\beta =\beta ^2`$ are, respectively, the mean photon numbers of the coupling and probe lasers in the coherent state $`|\alpha ,\beta `$. In the derivation of Eq. (11), we have used the large-$`n`$ approximation, i.e., $`\overline{n}_\alpha ,\overline{n}_\beta 1`$, so that the photon distributions are sharply peaked around their mean value. (For example, in a recent experiment of light storage , $`\overline{n}10^4`$.)

The polarization of the atomic medium at the frequency of the probe laser is determined to be

$$P(\omega _1)=\mu _{21}N\rho _{21}^A(\omega _1)=ϵ_0\chi (\omega _1)E_1(\omega _1),$$
(12)
where $`N`$ is the number density of atoms, $`\mu _{21}`$ the transition dipole moment between states $`|2`$ and $`|1`$, $`ϵ_0`$ the free space permittivity, $`\chi (\omega _1)`$ the susceptibility of the atomic medium, $`E_1(\omega _1)`$ the Fourier component of the mean electric field for the probe laser at frequency $`\omega _1`$. The value of $`E_1(\omega _1)`$ can be calculated from the ”steady” state (10). In the large-$`n`$ approximation, $`E_1(\omega _1)=_1\alpha `$, where $`_i=(\mathrm{}\omega _i/2ϵ_0V)^{1/2}`$ ($`i=1,2`$) with $`V`$ the quantized volume. Then the susceptibility is given by

$$\chi (\omega _1)=\frac{4N|\mu _{12}|^2\overline{\mathrm{\Omega }}_2^2\mathrm{\Delta }_1}{\mathrm{}ϵ_0(\overline{\mathrm{\Omega }}_1^2+\overline{\mathrm{\Omega }}_2^2)^2},$$
(13)
where $`\overline{\mathrm{\Omega }}_1=\mathrm{\Omega }_1(\overline{n}_\alpha ,\overline{n}_\beta )`$ and $`\overline{\mathrm{\Omega }}_2=\mathrm{\Omega }_2(\overline{n}_\alpha ,\overline{n}_\beta )`$ are the Rabi frequencies of the coupling and probe lasers, respectively. This is the main result of this paper, valid for arbitrary ratio of $`\overline{\mathrm{\Omega }}_1/\overline{\mathrm{\Omega }}_2`$.

Eq. (13) exhibits the signature of EIT: the linear susceptibility vanishes at the resonance ($`\mathrm{\Delta }_1=0`$). The derivative of $`\chi (\omega _1)`$ is related to the group velocity for the probe laser pulse: $`v_g=c/[1+(\omega _1/2)(d\chi /d\omega _1)]`$, with $`c`$ speed of light in vacuum. Thus we have

$$v_g=v_g^0\frac{\left(\overline{\mathrm{\Omega }}_1^2+\overline{\mathrm{\Omega }}_2^2\right)^2}{\overline{\mathrm{\Omega }}_2^4},$$
(14)
where $`v_g^0=\mathrm{}cϵ_0\overline{\mathrm{\Omega }}_2^2/(2\omega _1|\mu _{12}|^2N)`$ is the usual expression for the group velocity . This equation shows that in general the group velocity of the probe laser depends on the Rabi frequency of the coupling as well as the probe laser. The $`\overline{\mathrm{\Omega }}_1`$ dependence of the group velocity is a new result of our treatment, closely related to the higher order nonlinear susceptibilities we will discuss below.

To relate the present results to earlier studies of EIT-based dispersive properties , we expand the susceptibility given in Eq. (13) in terms of powers of $`\overline{\mathrm{\Omega }}_1/\overline{\mathrm{\Omega }}_2`$. When $`\overline{\mathrm{\Omega }}_1\overline{\mathrm{\Omega }}_2`$, the first term of this expansion gives us the linear susceptibility

$$\chi ^{(1)}(\omega _1)=\frac{4|\mu _{12}|^2N\mathrm{\Delta }_1}{\mathrm{}ϵ_0\overline{\mathrm{\Omega }}_2^2}.$$
(15)
Thus, at the lowest order we recover results in Ref. , when the decay rates of states are neglected at the resonant frequency of the probe laser. Similarly, keeping only the lowest-order term in $`\overline{\mathrm{\Omega }}_1/\overline{\mathrm{\Omega }}_2`$, Eq. (14) will give the usual expression, $`v_g^0`$, for the group velocity , which does not depend on $`\overline{\mathrm{\Omega }}_1`$ and involves the contributions only from linear susceptibility.

The refractive index of the medium can be obtained from Eq. (13) by definition $`n\sqrt{1+\chi }`$. Making use of Eqs. (13) and (14), we find the refractive index change near zero probe detuning to be

$$\mathrm{\Delta }n=\frac{\lambda _1}{2\pi }\frac{\mathrm{\Delta }_1}{v_g},$$
(16)
where $`\lambda _1`$ is the wavelength of probe laser. It is worthwhile to note that this refractive index change involves the contributions of all orders of nonlinear susceptibilities due to the intensity dependence of the group velocity in (14). For the slow light experiment , with $`v_g=17`$ m/s and parameters $`\mathrm{\Delta }_1=1.3\times 10^6`$ rad/s and $`\lambda _1=589`$ nm, we obtain from Eq.(16) $`\mathrm{\Delta }n=7.2\times 10^3`$. This value agrees with the measured value in Ref. .

Giant Non-linearities Nonlinearities play an important role not only in nonlinear optics but also in quantum optics. They may be used for generation of enhanced squeezing , quantum computation and quantum teleportation , and quantum nondemolition measurements . Since EIT takes place in the vicinity of atomic resonance, large nonlinearities are naturally expected. We note that in the conventional steady-state approach, it is difficult to obtain nonlinear susceptibilities higher than $`\chi ^{(1)}`$ and $`\chi ^{(3)}`$. However, in our formalism we can get arbitrary higher-order nonlinear susceptibilities once for all.

The higher-order nonlinear susceptibilities are defined by

$`\chi (\omega _1)`$ $`=`$ $`\chi ^{(1)}(\omega _1)+\chi ^{(3)}(\omega _1)|E(\omega _1)|^2`$ (18)
$`+\chi ^{(5)}(\omega _1)|E(\omega _1)|^4+\chi ^{(7)}(\omega _1)|E(\omega _1)|^6+\mathrm{}`$
where $`\chi ^{(1)}(\omega _1)`$ is the linear susceptibility given in Eq. (15), and $`\chi ^{(k)}(\omega _1)(k3)`$ represent the $`k`$th-order nonlinear susceptibility. Assuming that the coupling laser is stronger than the probe, we expand the susceptibility (13) in powers of $`\overline{\mathrm{\Omega }}_1/\overline{\mathrm{\Omega }}_2`$. Rewriting this expansion in the form of Eq. (18), we find

$`\chi ^{(3)}(\omega _1)`$ $`=`$ $`{\displaystyle \frac{4ϵ_0c}{I_2}}\chi ^{(1)}(\omega _1),\chi ^{(5)}(\omega _1)={\displaystyle \frac{3ϵ_0c}{I_2}}\chi ^{(3)}(\omega _1),`$ (19)
$`\chi ^{(7)}(\omega _1)`$ $`=`$ $`{\displaystyle \frac{8ϵ_0c}{3I_2}}\chi ^{(5)}(\omega _1),`$ (20)
where $`I_2=2ϵ_0_2c\beta ^2`$ is the intensity of the incident coupling laser.

Normally nonlinear optical properties of the medium are described by the nonlinear refractive indices :

$`n`$ $`=`$ $`n_0+n_2|E_1|^2+n_4|E_1|^4+n_6|E_1|^6+\mathrm{}`$ (21)
where the first nonlinear correction to the refractive index is the Kerr coefficient $`n_2`$, which is related to $`\chi ^{(3)}(\omega _1)`$. $`n_k(k4)`$ are higher-order nonlinear refractive-index coefficients, related to higher-order nonlinear susceptibilities up to $`\chi ^{(k+1)}(\omega _1)`$.

Making use of Eqs. (18) and (20), we obtain from Eq. (21) nonlinear refractive index coefficients:

$`n_2`$ $`=`$ $`{\displaystyle \frac{2ϵ_0c}{I_2}}\chi ^{(1)}(\omega _1)={\displaystyle \frac{2ϵ_0c\mathrm{\Delta }_1}{\pi I_2}}{\displaystyle \frac{\lambda _1}{v_g^0}},`$ (22)
$`n_4`$ $`=`$ $`{\displaystyle \frac{3ϵ_0cn_2}{I_2}},n_6={\displaystyle \frac{8ϵ_0cn_4}{3I_2}}.`$ (23)
To increase the value of nonlinear refractive-index coefficients one may either increase the atomic density or decrease the coupling laser intensity.

To demonstrate the magnitude of the giant nonlinearities derived above, we calculate the nonlinear refractive-index coefficients using the parameters ($`I_2=40`$ mW/c$`m^2`$, $`\mathrm{\Delta }_1=1.3\times 10^6`$ rad/s, and $`\lambda _1=589`$ nm) in the ultraslow light experiment reported in Ref. , in which a light pulse speed $`17`$ m/s was observed. We estimate that under these conditions, $`n_2=1.9\times 10^7`$ m2/V2, $`n_4=3.8\times 10^{12}`$ m4/V4, and $`n_6=6.7\times 10^{17}`$ m6/V6. In terms of a common practical unit, $`n_2=0.36\mathrm{cm}^2/\mathrm{W}`$, $`n_4=13.0\mathrm{cm}^4/\mathrm{W}^2`$, and $`n_6=4.5\times 10^2\mathrm{cm}^6/\mathrm{W}^3`$. The value of the Kerr nonlinearity is of the same order of magnitude as that indirectly measured in Ref. , almost $`10^6`$ times greater than that measured in cold Cs atoms , and $`10^{12}`$ times greater than that measured in other materials . The fourth-order refractive-index coefficient $`n_4`$ is $`10^{22}`$ times greater than that measured in other materials .

In optical fibers, the ratio between the second- and the fourth-order refractive-index coefficients is an essential parameter to obtain stable spatial solitary waves. The lower the ratio $`n_2/n_4`$, the lower the required power for stable beam propagation. For the atomic medium with EIT, we have $`n_2/n_4=I_2/(3ϵ_0c)`$. With the parameters in the experiment, we estimate $`|n_2/n_4|10^2\mathrm{W}/\mathrm{cm}^2`$. This ratio is small compared to most other nonlinear media  by almost 11 orders of magnitude. From Eq. (23) we also obtain the ratio between the fourth- and the sixth-order refractive-index coefficients $`n_4/n_6=3I_2/8(ϵ_0c)`$, which is of the same order as $`n_2/n_4`$.

We note that when the ratio of $`\overline{\mathrm{\Omega }}_1/\overline{\mathrm{\Omega }}_2`$ is close to unity, one had better directly deal with the original formula (10), rather than using its power series expansion.

Conclusions We have developed a fully quantum treatment of EIT in atomic medium. Both the probe and coupling lasers are quantized and treated on the same footing. This allows us to deal with the cases in which the ratio of the probe-to-coupling Rabi frequencies is not necessarily small. At the lowest order in this ratio, we are able to reproduce the known results for slow pulse propagation. At higher orders we have uncovered that accompanying the EIT phenomenon, atomic medium possesses giant optical non-linearities, which can give rise to dramatic enhancement of the Kerr as well as higher order refractive-index coefficients. In other words, the atomic medium with EIT is really a very unusual optical medium; to the list of its astonishing optical properties, we now add another remarkable one: giant optical non-linearities. It would be very interesting to observe these giant non-linearities and explore their potential applications.

Contrary to the usual treatment of EIT, our present treatment is essentially a time-independent approach, because both probe and coupling lasers are included as part of the dynamical system. On the other hand, an important limitation is that we have ignored the decay rates of various levels. How to incorporate them in a fully quantum treatment is a challenging issue.

This work was supported in part by the US NSF under grant PHY-9970701, and a Seed Grant from the University of Utah. L.M.K. also acknowledges support from the Climbing Project, NSF and EDEYTF of China, and Hunan Province STF.