A Bartecka
email: bartecka@uni.opole.pl

A Bac Lawski

J Musielok

Stark broadening studies of the infrared OI multiplet 3p 5 P -3d 5 D o

Keywords: numbers: 32, 70, Jz, 52, 70, Kz

Experimental line broadening studies of the infrared OI multiplet 3p 5 P -3d 5 D o are reported. Stark broadening parameters: electron impact widths (w) and shifts (d) as well as asymmetry (A) are determined applying a wall-stabilized high-current arc as the plasma source. Broadening studies are performed in the electron density range 0.8•10 16 ÷1.4•10 16 cm -3 , corresponding to plasma temperatures from 12000 to 13800 K. The obtained broadening and shift parameters, normalized to the standard electron density value of 1•10 16 cm -3 , are compared with results of calculations reported in the monograph by Griem [Spectral Line Broadening Plasmas,

Introduction

Line shape and shift studies of spectral transitions of neutral oxygen (OI) have been the subject of numerous experimental and theoretical papers. Critical reviews of experimental data on OI Stark widths and shifts, published in the time period from the early sixties up to the end of the twentieth century, are reported by: Konjević and Roberts [1], Konjević and Wiese [2], and Konjević at al. [3]. The comparison of experimental data with Stark broadening and shift parameters calculated by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF] shows a differentiated picture: for some multiplets the agreement is satisfying, but in other cases the discrepancies are clearly outside the uncertainties of measurements. For example, in the case of prominent visible (red) quintet transition 3p 5 P -4d 5 D o at wavelength 615.73 nm, the old study based on wall-stabilized arc emission measurements by Jung [START_REF] Jung | [END_REF], yields line widths about 3 times smaller than those calculated by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF]. On the other hand, gas-driven shock tube measurements by Miller and Bengtson [6] provide Stark widths somewhat larger (20%) than those resulting from calculations [START_REF] Griem | Spectral line broadening by plasmas[END_REF]. In the case of higher excited multiplet 3p 5 P -5d 5 D o at 532.9 nm, the two so far reported measurements provide narrower line widths than the calculations of Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF] -the shock tube results [6] by about 20%, while the arc emission experiment of Goly and Weniger [7] even by 40%. More recent calculations of Stark broadening parameters performed by Ben Nassib et al. [8] and Rivièrie [9] provide also, among others, electron impact line widths for the OI multiplet at 529.9 nm. The newly calculated data for this multiplet are in good mutual agreement and are indeed by about 20% narrower than those given by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF]. However, the calculated widths are still by about 30% larger than the experimental data reported by Goly and Weniger [7]. Unfortunately, both recent calculations do not provide data for the less excited multiplets: 3p 5 P -4d 5 D o at 615.7 nm and 3p 5 P -3d 5 D o at 926.4 nm. In the monograph by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF], Stark broadening data are reported also for the transition with the longest wavelength in the series of OI quintet transitions 3p 5 Pnd 5 D o , namely for the transition 3p -3d. The electron impact width calculated by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF] for this transition is significantly narrower than for the above discussed higher excited multiplets. This broadening parameter for the 3p -3d line is about 5 and 13 times narrower than for the 3p -4d, and 3p -5d transitions, respectively. According to our knowledge, experimental studies of Stark broadening for the infrared multiplet have not been performed so far. Recent development in the construction of charge coupled devices made the registration of spectra in the infrared part of the spectrum considerably easier. Therefore, we decided to undertake Stark broadening studies for this quintet infrared transition at wavelength 926.4 nm.

Preliminary analysis of the multiplet 3p

5 P -3d 5 D o
The quintet transition 3p 5 P -3d 5 D o consists of nine spectral lines spread in the spectrum from 926.0806 nm up to 926.6006 nm. In table 1 we present the spectroscopic data of this multiplet taken from the NIST data base [START_REF]NIST Atomic Spectra Database[END_REF]. Quantum transitions with the same lower level appear very close in wavelength -the components are spread by 13.0, 19.4 and 18.0 pm, for quantum transitions ending up on the levels J = 1, 2 and 3, respectively. In the case of line broadening caused by charged particles acting on the emitter in plasmas, these individual components, having the same final lower level, usually cannot be resolved. Therefore, in the case of spectra emitted from plasmas, normally only three spectral features are observed with peaks at 926.087 nm, 926.272 nm and 926.599 nm. In order to facilitate further studies of the measured spectra we performed a preliminary analysis of the three above mentioned peaks, formed by closely spaced individual components, broadened by the Stark effect. Stark profiles have been calculated applying the well known j(x) formula given by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF] and assuming different broadening parameters: electron impact widths w and asymmetry parameters A, corresponding to different electron densities in the plasma. Assuming that each spectral line (fine structure component) has the same profile, the shapes of each peak have been calculated taking into account the unperturbed wavelength positions of each fine structure component and their relative intensities obtained from line strengths listed in column 2 of table 1. The shapes of the three peaks obtained in this way have been then fitted with j(x) profiles and the resulting ("joint") broadening parameters w and A have been determined. As an example, in figure 1 we show the three fine structure components forming the strongest peak of the multiplet at 926.599 nm (dashed lines), their sum (circles) and the fitted Griem-like profile for two pairs of "input" data: w = 5 pm, A = 0.05 (upper plot) and w = 40 pm, A = 0.2 (lower plot). The best fits (solid lines) correspond to: w = 6.7 pm, A = 0.00013 (upper plot) and w = 40.4 pm, A = 0.198 (lower plot). As can be seen from the upper graph of figure 1, the wavelength positions of the two weaker components enhance the violet wing of the peak causing an artificial asymmetry. Thus, the fit of the Griem-like profile to the sum (peak) seems not to be perfect. The determined asymmetry parameter A is obviously substantially too small and also the resulting electron impact width is significantly overestimated. At higher electron densities of the plasma (larger w values), the artificial asymmetry caused by the two weak components on the blue wing does not significantly influence the fitting procedure. As can be seen from the lower plot of figure 1, the Griem-like profile satisfactorily describes the "true" asymmetry resulting in the enhanced red wing of the peak and the "joint" data w and A do not differ substantially form the corresponding "input" data. The results of this preliminary analysis for all three peaks are shown in figure 2. In the upper part of this figure the ratios w joint /w input are shown as function of w input , in the lower part the analogous graph is presented for the asymmetry parameter. As can be seen, in the case of two peaks (with the shortest and the longest wavelength) for "input" values w ≥ 25 pm the "joint" values for the peaks are only somewhat larger (less than 4%) than the input data for individual components. In the case of the central peak at λ = 926.272 nm, the discrepancy is systematically about two times larger. The reason for that is simple -the wavelength difference between the two strongest components exceeds 10 pm and the wavelength interval between the outermost components amounts 19.4 pm, leading to significantly larger "observed" widths of this peak. By treating three line components as an individual spectral feature (peak), the true asymmetry of a component is obviously somewhat disturbed by the wavelength structure of this peak. Indeed, as can be seen from the lower part of figure 2, the "joint" (observed) asymmetry parameter is systematically smaller than the corresponding input value, with the exception of only very low A "input" values in the case of the peak with the shortest wavelength. However, for input A values exceeding 0.15, all three peaks exhibit "joint" (observed) asymmetries which are systematically by about 10% smaller than the input values. In this analysis, the influence of the Doppler and instrumental broadening has not been taken into account. As mentioned in section 3, the instrumental profile could be well approximated by a Gaussian profile with a Full Widths at Half Maximum (FWHM) of 8 pm. The FWHM resulting from the Doppler effect amounts between 17 and 19 pm, depending on the plasma temperature. Both effects, the instrumental and the Doppler broadening, lead to a final Gaussian with a FWHM of 20 to 22 pm. The wavelength differences between the fine structure components forming a given peak are of the same order as the resulting convolution of the Doppler and instrumental widths. Therefore, in case of including this "additional" broadening, the calculated ratios w joint /w input and A joint /A input shown in figure 2 are expected to be systematically somewhat closer to the desired value of one. Indeed, when the Stark profile of each fine structure component was convoluted with the convolution of the Doppler and instrumental profile, the obtained ratios w joint /w input and A joint /A input were systematically slightly closer to one, particularly for smaller values of the "input" data of w and A.

Experimental procedure

A wall-stabilized high-current arc operated at atmospheric pressure was applied as the plasma source. The scheme of the arc construction, including the scheme of the gas injection system, is shown in figure 3. The arc column (85 mm long, diameter 4 mm) was operated in a mixture of helium (gas flow rate 45 cm 3 /s), argon (3.4 cm 3 /s), carbon-dioxide (0.78 cm 3 /s) and hydrogen (0.51 cm 3 /s). The arc regions in the direct vicinity of both carbon electrodes were additionally fed with argon at a flow rate of 2.1 cm 3 /s. This rinsing argon stream protects the carbon electrodes from erosion. The amount of CO 2 and hydrogen was carefully adjusted in order to obtain optimal conditions for registration of (i) the hydrogen H β line (for electron density determination) and (ii) selected oxygen spectral lines (the studied multiplet at 926 nm and a set of "thermometric" transitions). At the chosen gas flow rates and arc currents of 39 and 54 A, plasma conditions (electron densities and temperatures) were achieved allowing accurate Stark broadening measurements of the OI multiplet 3p 5 P -3d 5 D o to be performed. At these conditions the three main spectral peaks do not substantially overlap and thus, Stark broadening studies were possible.

In figure 4 the scheme of the experimental setup is shown. The arc (1) was imaged in end on direction (magnification 1:1) via the swivel plane mirror ( 5) and the concave mirror (4) (focal length f = 730 mm) on the entrance slit (26 µm) of the Ebert-type grating spectrometer ( 6) equipped with a CCD array detector Symphony (7) (pixels 26×26 µm) and the computer (8). By changing the angular position of the mirror (5) by 180 o the radiation originating from the tungsten strip standard source (2) was registered in all studied spectral intervals. In this way, the measured intensities originating form the arc source could be converted to absolute radiance units. The concave mirror (3) was applied for performing self-absorption checks. By reflecting the radiation reaching this mirror back into the arc and comparing the measured signals, with and without back reflection, conclusions about possible self-absorption effects in the arc plasma could be drawn. For our plasma conditions, these checks demonstrated negligible self-absorption in all studied spectral intervals (OI multiplets and the H β transition). The spatial resolution of the optical imaging system perpendicular to the arc axis, allows selecting the radiation originating from plasma layers parallel to the arc axis of a size of 0.1 mm. The CCD detector area was divided into 64 tracks, each consisting of four pixels in the direction perpendicular to the arc axis and 1024 pixels along the direction of the spectrometer's dispersion. The distance between the centers of neighboring tracks amounts 0.1 mm and precisely matches the spatial resolution of the optical imaging system. In this manner, each detected spectrum yields directly the radial intensity distribution in the arc plasma for the selected spectral interval. The spectrometer (6) was used in the double path mode applying a grating with 1300 grooves/mm. This instrumentation allows to obtain spectra with reciprocal dispersion of 4.3 pm/pixel and 4.8 pm/pixel, depending on the studied wavelength range. Other details, concerning the arc plasma source, the optical imaging system and the spectroscopic instrumentation, can also be found e.g. in [START_REF] Wujec | [END_REF]12,13].

Plasma diagnostics

As mentioned above, two experiments were performed at arc currents of 39 and 54 A. For each arc current, the radiation emerging in end on direction was registered in four different spectral intervals around the following wavelengths: 486.1 nm (hydrogen H β line), 777.3, 882.0 and 926.4 nm (selected OI transitions). Each spectral interval was registered at least three times and the results were averaged. From measured OI spectra the total line (multiplet) intensities have been obtained by integrating the respective spectral intervals and calibrating them against the respective intensities of the standard source. The quintet transitions at 777.3 and 926.4 nm originate from upper levels with excitation energies of 10.741 and 12.079 eV, respectively, while the singlet line at 882.0 nm arises from the upper level 3p 1 F 3 with E k = 14.134 eV. Using the corresponding transition probabilities (taken from the NIST data base [START_REF]NIST Atomic Spectra Database[END_REF]) and applying the Boltzmann plot method, the radial distributions of the arc temperature for the two experiments at 39 and 54 A have been determined. The results are presented in the upper part of figure 5.

The main plasma parameter determining the Stark broadening of spectral lines is the electron density. This crucial plasma parameter was determined from measured FWHM of the hydrogen H β transition. For this purpose, we applied the theoretical broadening results of Gigosos and Cardeñoso [14]. These calculations include the contribution to the observed line widths arising from ion dynamic effects. Other broadening mechanisms, such as instrumental and Doppler broadening, have been taken into account while determining the pure Stark FWHM of the H β line.

The low pressure Plücker tube run in hydrogen was used as a light source for determination of the instrumental profile in this spectral region. It turns out that the instrumental profile could be well approximated by a Gaussian with a FWHM of 8 pm. The contributions to the H β line widths originating from the Doppler broadening amounts from 34 to 40 pm, depending somewhat on the plasma temperature. Thus, the total contribution from both effects (instrumental and Doppler broadening) amounts from 35 to 41 pm. These widths are significantly smaller (at least by a factor of 40) than the measured FWHM of the H β line. Nevertheless, evaluating the pure Stark FWHM the necessary corrections of directly measured FWHM of the H β line have been performed. Using the obtained pure Stark widths and the theoretical broadening data for the H β line, published by Gigosos and Cardeñoso [14], the electron density distributions have been evaluated. The results are shown in the lower part of figure 5. As can be seen, in the central part of the arc (r≤0.9 mm and r≤0.5 mm, for experiments at 54 A, and 39 A, respectively), the radial gradients of both main plasma parameters are small. In the case of both applied arc currents, the gradients of the temperature (∂T /∂r) do not exceed the value of -1200 K/mm for these plasma volumes. The corresponding gradients of the electron density, for the two arc currents 54 and 39 A, do not exceed:

(∂n e /∂r) = -7•10 15 cm -3 /mm and (∂n e /∂r) = -2.6•10 15 cm -3 /mm, respectively. As a result of these small gradients, within the distance of 0.1 mm in the arc (corresponding to one track) the temperature (T ) and electron density values (n e ) differ by not more than 0.5% and 4% of the respective mean T and n e values. Thus, the spectra measured from these parts of the arc can be regarded as originating from nearly homogeneous plasma layers.

Results and discussion

As mentioned earlier, the wavelength interval embracing the quintet OI multiplet 3p 5 P -3d 5 D o was registered three times and the results were averaged. The long term stability of the arc as well as the reproducibility of the measured individual spectra was satisfactory. The differences between the parameters of individual spectra (intensities and widths) do not exceed 3%". To determine the Stark broadening parameters of this multiplet only the spectra originating from nearly homogeneous plasma layers have been applied, i.e. emitted from arc volumes with radii not exceeding 0.5 mm and 0.9 mm, for experiments at currents 39 and 54 A, respectively. For these selected plasma layers the observed "joint" widths (half widths at half maximum -HWHM) of the three peaks exceed the value of w = 25 pm and A = 0.15. According to the results of preliminary analysis shown in figure 2, it is possible to make reliable corrections of measured (observed) peak widths and asymmetry parameters to determine the respective quantities of the single fine structure component. Each measured spectrum (after radiance and wavelength calibration), taken from given plasma layer of known electron density and temperature, was fitted according to the following procedure:

1. The Griem-like profile (converted to wavelength units) was generated according to the formula [15]:

j A,R (λ) = j max ∞ 0 H R (β) 1 + [((λ -λ 0 )/w) -A 4/3 β 2 ] 2 dβ (1)
where: j max is the peak maximum, w is the electron impact broadening parameter (HWHM), A is the asymmetry parameter, H R (β) is the ion microfield distribution function at the (neutral) emitter position taken from [16]. R = r 0 /D is the so-called screening parameter (r 0 = (4πn e /3) -1/3 , D -Debye radius), which is needed for obtaining the proper microfield distribution function, β = F/F 0 is the reduced field strength (F 0 is the normal field strength depending on electron density). 2. This profile was numerically convoluted with a Gaussian being the outcome of the convolution of the instrumental and Doppler broadening.

3. A background continuum in the form C = a • λ + b was added to the three Griem-like profiles convoluted with the Gaussians and fitted to each measured spectrum of the multiplet.

The fitting procedure anticipates 10 free parameters to be determined: 3 wavelength positions (λ 0 ) and 3 maximum intensity values (j max ) for the three peaks, the parameters a and b, determining the course of the background, the electron impact broadening parameter w and the asymmetry parameter A for the multiplet, i.e., the same w and A values for all three peaks. In this procedure, the sum of squared differences between the measured and the calculated spectrum, taken over the whole wavelength interval was minimized. In figure 6, such spectra corresponding to two selected plasma layers, taken from experiments 39 and 54 A are shown as examples.

Within our fitting procedure we assumed that the shapes of all three peaks are identical (have the same w and A). Therefore the fitted peak intensities (j max ) are proportional to total intensities of the peaks and may be compared with relative line strengths taken from the NIST data base [START_REF]NIST Atomic Spectra Database[END_REF]. As mentioned earlier, each peak is formed by three fine structure components contributing to the observed peak intensity. The line strength of a given peak is simply equal to the sum of the strengths of the contributing fine structure components. Thus, the strengths of successive peaks (with rising wavelengths) obey the following ratios 1.00 : 1.67 : 2.33. Since the measurements yield intensity ratios, not line (peak) strength ratios, two corrections are required arising from: (i) the Boltzmann factor exp (E 2 -E 1 )/kT , and (ii) the fourth power of wavelength ratios (λ 2 /λ 1 ) 4 , in order to compare the theoretical and experimental data. Since the excitation energies and wavelengths are nearly the same, both these corrections are negligibly small -the Boltzmann factor, and the "wavelength" factor yield corrections not exceeding 0.1%. The ratios of the fitted peak values j max perfectly agree with those resulting from the LS coupling listed in table 1 -the discrepancies do not exceed 0.2%.

Shift measurements (d)

A glance on the spectra shown in figure 6 reveals a clear red shift of the multiplet, particularly at the higher electron density condition. However, this apparent observation could be quantitatively more or less uncertain. The second experiment at 54 A was performed after turning back the spectrometer grating to the appropriate position (wavelength range). Therefore, the wavelength calibrations of the spectra, performed against the NeI spectral line at 927.55191 nm, for both arc currents have been performed independently. Thus, at the apparatus width of our spectral instrumentation (FWHM = 8 pm, about 2 pixels), uncertainties of the wavelength calibration of the same order have to be expected. Moreover, as discussed in section 2, the wavelength structures of individual peaks introduce additional uncertainties in determination of the multiplet shift as a whole. Nevertheless, despite possible large errors of our shift measurements we decided to include this data in this paper. As mentioned above, the fitting procedure yields e.g. the wavelength positions of the three observed peaks. In order to obtain the shifts of the three peaks, their unperturbed wavelength positions (in air) have been calculated. Since each peak consists of three fine structure components, we defined the unperturbed peak wavelengths as:

λ peak = Σ k (λ ki •S ik )/Σ k (S ik ).
The wavelengths (λ ki ) and the line strengths (S ik ) listed in column 2 of table 1 have been taken for evaluating the peak wavelengths. In this way, the following unperturbed peak wavelengths have been obtained: 926.0870, 926.2718 and 926.5986 nm. The shifts for each peak have been calculated using the corresponding unperturbed wavelengths as reference values. At given plasma conditions, the determined shifts for individual peaks are nearly the same (the differences do not exceed the wavelength corresponding to 1 2 of a pixel), therefore we decided to show only the averaged values, which have to be considered as shift values of the multiplet as a whole.

In the upper part of figure 7 the shifts determined in this way are shown as a function of electron density of the plasma, while in the lower part, the shifts normalized to the standard electron density value of 10 16 cm -3 , are plotted against the temperature of the plasma. As mentioned earlier, the uncertainties of our shift measurements in the absolute scale are estimated to reach even 8 pm, i.e., of the same order as the measured shift values. Nevertheless, our measurements confirm the expected linear increase of the shift with increasing electron density. The averaged shift for the multiplet, normalized to the standard electron density of 10 16 cm -3 , amounts about 10 pm and is about 2 times smaller than the shift calculated by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF].

Determination of the electron impact broadening parameter (w)

In the upper part of figure 8, we plot the determined half widths at half maximum (HWHM) in pm (w) for the three peaks as a function of electron density. In the lower part of this figure the HWHM, reduced to the standard electron density value of 10 16 cm -3 (w nor ), are shown as a function of the plasma temperature. As can be seen from the upper graph of figure 8, the measured electron impact widths (w) -as expectedincrease linearly with increasing electron density. The reduced widths (w nor ), shown in the lower graph, slightly decrease with rising temperature. Our analysis described in section 2 leads to the conclusion that our measured peak widths are slightly larger than the widths of the fine structure components. This overestimation somewhat depends on the widths itself -larger corrections are required for smaller w values. Applying the results of this analysis, we corrected our measured peak widths in order to obtain the widths of the fine structure components. The dotted-dashed line represents the electron impact widths (HWHM) for the fine structure components. On average, our normalized electron impact width is 28±5 pm and is about 20% larger than the calculated result by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF]. In accordance with theoretical expectations, our measurements confirm the linear dependence of both Stark broadening parameters -the electron impact width (w) and the shift (d) on electron density. In figure 9, the ratios of these two quantities are shown versus the electron density. Indeed, the ratios d/w are the same within the uncertainty limits of our experiment in the range of the electron density from 0.8•10 16 to 1.4•10 16 cm -3 . The averaged ratio d/w amounts 0.38 and is about 2.5 times smaller than the result of Griem's calculations.

Studies of the asymmetry parameter (A)

In figure 10 we show our results determined for the dimensionless asymmetry parameter (A). Similarly as in figure 8, in the upper part the asymmetry parameter is plotted against the electron density, while in the lower part, the reduced asymmetry parameter versus the plasma temperature. Again a weak linear dependence of the asymmetry parameter on electron density is found. The reduced A nor values slightly increase with rising temperature. In case of the asymmetry parameter the correction discussed in section 2 leads to somewhat larger A values for the fine structure components. As can be seen from figure 10 (dotted-dashed line), the correction does not exceed 10%. The comparison with the result reported by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF] also reveals large discrepancy -our reduced asymmetry parameter is about 2.5 times larger.

Concluding remarks

The determination of Stark broadening parameters for the fine structure components of the multiplet 3p 5 P -3d 5 D o is significantly hindered by the complicated structure of this transition. Three sets of fine structure components, each consisting of three transitions, form three rather good isolated peaks at lower electron densities. At higher electron densities these peaks overlap substantially. Therefore, the determination of Stark broadening parameters for the individual fine structure component is a difficult task. Among the three Stark broadening parameters studied in this work: the electron impact shift (d), the electron impact width (w) and the asymmetry parameter (A), the most reliable results are obtained for the electron impact width. This parameter agrees -within the error limits of our experiment (20%) and the uncertainty of calculations (20%) -with the result reported by Griem [START_REF] Griem | Spectral line broadening by plasmas[END_REF]. In case of the two other parameters d and A, the agreement between measured and calculated data is less satisfying. Large discrepancy, exceeding a factor of 2, is found in the case of the asymmetry parameter. The uncertainties of our shift measurements in the absolute scale do not allow to make indisputable conclusions about the grade of consistence between experiment and theory. Nevertheless, our "relative" shift measurements, i.e., the quantitative increase of the shift with increasing electron density are trustworthy.

Figure captions

Figure 1 .

 1 Figure 1. The three fine structure components forming the strongest peak of the quintet multiplet at 926.599 nm (dashed lines), their sums (circles) and the fitted Griem-like profiles for two pairs of "input" data are shown: w = 5.0 pm, A = 0.05 (upper plot) and w = 40.0 pm, A = 0.20 (lower plot).

Figure 2 . 5 Figure 3 .

 253 Figure 2.Electron impact broadening (w) and asymmetry (A) parameters of individual fine structure components ("input") and of the observed spectral peaks ("joint") are compared. Upper graph -the ratio w joint /w input is shown a function of width (in pm) of the individual component w input , and lower graph -the ratio A joint /A input is shown versus the dimensionless asymmetry parameter A of the individual component A input . The resulting "joint" profiles were fitted using the j(x) functions with broadening parameters w joint and A joint (both free) -for details see the text in the section 2

9 Figure 4 .

 94 Figure 4. The scheme of the experimental setup is shown: 1 -wall-stabilized arc, 2radiometric standard source, 3 -and 4 -concave mirrors, 5 -swivel plane mirror, 6 -Ebert-type grating spectrometer, 7 -CCD detector, 8 -computer, 9 -spectral filter.

Figure 5 .

 5 Figure 5. The radial distributions of the temperature (upper graph) and electron density of the plasma (lower graph) are shown. Stark broadening studies of the OI multiplet have been performed on the basis of spectra obtained from arc volumes with radii less than 0.9 mm (expt. 59 A) and 0.5 mm (expt. 39 A), as indicated by the arrows.

Figure 6 .

 6 Figure6. The spectral intensity distributions within the OI multiplet 3p 5 P -3d 5 D o corresponding to two plasma conditions are shown. The squares represent the measured intensities after radiance and wavelength calibration, the lines -the results of the fitting procedure. The three short vertical lines mark the positions of the unperturbed wavelengths of the three peaks. These positions are "weighed" wavelengths of the corresponding fine structure components (for details see text).

Figure 7 .

 7 Figure 7. The measured absolute red shifts (d) of the multiplet as a function of the electron density (upper graph), and the shifts normalized to the standard electron density of 10 16 cm -3 (d nor ) as a function of the plasma temperature (lower graph), are shown. The error bars shown in both graphs do not include the possible systematical error originating from the wavelength calibration in an absolute scale. This systematical error may reach even 8 pm.

Figure 8 .

 8 Figure 8. Upper graph -the determined half widths at half maximum (w) in pm for the three peaks are plotted as a function of electron density; lower graph -the reduced (to the standard electron density value of 10 16 cm -3 ) HWHM (w nor ) (in pm) are presented as a function of the plasma temperature. The dashed lines represent the best fit to the experimental points, the dotted-dashed line in the lower graph represents the best fit results after corrections resulting from the analysis presented in figure2,i.e., the final w nor values, corresponding to the individual fine structure component of the multiplet. The dotted line in the lower graph represents the theoretical data reported by Griem[START_REF] Griem | Spectral line broadening by plasmas[END_REF].

Figure 9 .

 9 Figure 9. The measured ratios: electron impact shift/electron impact width, (d/w), are shown versus the electron density of the plasma. As expected, this ratio not depend on the electron density. The averaged value is 0.38.

Figure 10 .

 10 Figure10. Upper graph -the determined dimensionless asymmetry parameters are shown as a function of the electron density; lower graph -the reduced asymmetry parameters (to the standard electron density value of 10 16 cm -3 ) (A nor ) are plotted versus the plasma temperature. The dashed lines represent the best fit to the experimental points, the dotted-dashed line in the lower graph represents the best fit results after corrections resulting from the analysis presented in figure2, i.e., the final A nor values corresponding to the individual fine structure component. The asymmetry parameter A nor , calculated by Griem[START_REF] Griem | Spectral line broadening by plasmas[END_REF] (dotted line), slightly decreases with temperature and is -on average -about 2.5 times smaller than our data.

Table 1 .

 1 Selected spectroscopic data concerning the OI multiplet 3p 5 P -3d 5 D o λ(nm) S ki (a.u.) E i (cm -1 ) E k (cm -1 ) J i J k

	926.0806 17.5	86625.76	97420.99	1	0
	926.0848 39.3	86625.76	97420.94	1	1
	926.0936 30.6	86625.76	97420.84	1	2
	926.2582 13.1	86627.78	97420.94	2	1
	926.2670 51.0	86627.78	97420.84	2	2
	926.2776 81.6	86627.78	97420.72	2	3
	926.5826 5.84	86631.45	97420.84	3	2
	926.5932 40.7	86631.45	97420.72	3	3
	926.6006 157	86631.45	97420.63	3	4


P Jönsson

T Carette

M Nemouchi

M Godefroid

Ab initio calculations of 14 N and 15 N hyperfine structures

Keywords: numbers: 31.15.ac, 31.30.Gs, 32.10.Fn Electron correlation, hyperfine structures, Nitrogen spectra

Introduction

Doppler-free spectra of the near-infrared N I transitions in the 2p 2 ( 3 P )3s 4 P → 2p 2 ( 3 P )3p 4 P o and 2p 2 ( 3 P )3s 4 P → 2p 2 ( 3 P )3p 4 D o multiplets have been recorded by Jennerich et al. [1] using saturated absorption spectroscopy, extending the pioneer work of Cangiano et al. [2] using a similar set-up. From the analysis of these spectra, Jennerich et al. [1] determined the hyperfine structure constants for the various J-values of the three multiplets involved, for both isotopes 14 N and 15 N. The isotope shifts in each multiplet have also been measured, revealing a significant J-dependence of the shifts. These authors recommended a theoretical investigation of the underlying cause of this unexpected phenomenon. Similar measurements, using Doppler-free saturated absorption spectroscopy, have been performed in atomic Fluorine [3], Chlorine [4] and Oxygen [5].

The present work presents a robust ab initio theoretical estimation of the relevant hyperfine structure parameters, using the ATSP2K package [6], based on the non relativistic multiconfiguration Hartree-Fock (MCHF) and configuration interaction (CI) methods. The minimum theoretical background is presented in section 2 for the hyperfine interaction. Section 3 is dedicated to the description of the electron correlation models. A refined calculation for the even-parity term is presented in section 4. A first comparison between theory and observation appears in section 5, revealing unexpected discrepancies. Relativistic effects are investigated in section 6 through the relativistic configuration interaction (RCI) approach, confirming serious problems in the experimental analysis as discussed in section 7.

The hyperfine interaction

The theory underlying MCHF calculations of hyperfine structure parameters can be found in [7,8]. Neglecting the relativistic effects, the diagonal and off-diagonal A and B hyperfine interaction constants are expressed in terms of the J-independent orbital (a l ), spin-dipole (a sd ), contact (a c ) and electric quadrupole (b q ) electronic hyperfine parameters defined by

a l = γLSM L M S | N i=1 l (1) 0 (i)r -3 i |γLSM L M S , (1)
a sd = γLSM L M S | N i=1 2C (2) 0 (i)s (1) 0 (i)r -3 i |γLSM L M S , (2)
a c = γLSM L M S | N i=1 2s (1) 0 (i)r -2 i δ(r i )|γLSM L M S , (3)
b q = γLSM L M S | N i=1 2C (2) 0 (i)r -3 i |γLSM L M S , (4)
and calculated for the magnetic component M L = L and M S = S. The first three parameters (1), ( 2) and (3), contribute to the magnetic dipole hyperfine interaction constant through

A J = A orb J + A sd J + A c J , (5)
with

A orb J = G µ µ I I a l L. J LJ(J + 1) , (6)
A sd J = 1 2 G µ g s µ I I a sd 3 L. S L. J -L(L + 1) S. J SL(2L -1)J(J + 1) , (7)
A c J = 1 6 G µ g s µ I I a c S. J SJ(J + 1) , (8)
while the last one (b q ) constitutes the electronic contribution to the electric quadrupole hyperfine interaction

B J = -G q Q b q 6 L. J 2 -3 L. J -2L(L + 1)J(J + 1) L(2L -1)(J + 1)(2J + 3) . ( 9
)
Expressing the electronic parameters a l , a sd and a c in atomic units (a -3 0 ) and µ I in nuclear magnetons (µ N ), the magnetic dipole hyperfine structure constant A is calculated in units of frequency (MHz) by using G µ = 95.41067. Similarly, the electric quadrupole hyperfine structure constant B is expressed in MHz when adopting atomic units (a -3 0 ) for b q , barns for Q and G q = 234.96475. The electronic parameter governing the mass isotope shift of an atomic energy level is the S sms parameter [9] defined by

S sms = -γLSM L M S | N i<j ∇ i • ∇ j | γLSM L M S . (10)

Correlation models

The multiconfiguration Hartree-Fock (MCHF) variational approach consists in optimizing the one-electron functions spanning a configuration space and the mixing coefficients of the interacting configuration state functions (CSF) [START_REF] Fischer | Computational Atomic Structure -An MCHF Approach[END_REF] for describing a given term

Ψ(γLSM L M S ) = i c i Φ(γ i LSM L M S ). ( 11
)
Efficient MCHF expansions are often built by allowing single and double excitations from a multireference set (SD-MR-MCHF). As far as hyperfine structures are concerned, successful applications of this method are found for light elements such as Li-like ions [9], Be I, B I, C II and C I [START_REF] Jönsson | [END_REF]12], N I [13], O I [13,14] or Na I [15].

For the even parity term 2p 2 3s 4 P , the configuration expansion is generated from single and double (SD) excitations from the {2s 2 2p 2 3s, 2p 4 3s, 2s2p 4 } (12) multireference (MR) to increasing active sets (AS) of orbitals that are denoted by specifying the number of orbitals for each l-symmetry. This multireference set captures the dominant correlation effects through a physical good "zero-order" wave function. We include only configuration state functions (Φ(γ i LSM L M S )) that interact with the multireference, adopting the reversed orbital order, i.e. coupling sequentially the subshells by decreasing n and l. This technique reduces substancially the size of the MCHF expansions while keeping the dominant correlation contributions [12,16]. With the largest set of orbitals (10s9p8d6f 3g) optimized through these calculations, the effect of higher excitations is investigated through configuration interaction (CI) calculations using the configuration state function set obtained by adding to the original SD-MR CSF list, the triple and quadruple excitations from the same multireference to smaller orbital active sets (up to 6s5p4d3f ). For generating these lists, some limitations have been introduced by imposing the restriction that there should be at least 5 orbitals with n ≤ 3 in the CSFs produced. The merging of the original SD-MR and the TQ-MR CSF lists is noted by the union (∪) symbol.

For the odd parity states 2p 2 3p 4 P o and 2p 2 3p 4 D o , a similar strategy is adopted, using the following multireference set

{2s 2 2p 2 3p, 2p 4 3p, 2s2p 2 3p3d} . ( 13
)
The total energies (E), the specific mass shift (S sms ) and the hyperfine interactions parameters {a l , a sd , a c , b q } defined in section 2, are reported in Tables 1, 2 and 3, respectively. These tables illustrate the smooth convergence of the various parameters while improving the correlation model by increasing the size of the orbital active set. For the even parity term 2p 2 3s 4 P (Table 1), triple and quadruple excitations from the multireference affect the a l and a sd parameters at the level of 2% while the contact term a c is much more sensitive (33%). For the odd parity 2p 2 3p 4 P o and 2p 2 3p 4 D o terms (Tables 2 and3), all hyperfine parameters reach a high degree of convergence, except the contact parameter that is strongly affected (around 30%) by TQ excitations. For the three terms, the last layer added at the SDTQ-CI level of approximation brings a negligible increment, except for the contact contributions (3%). The convergence patterns of the contact contributions, as well as the sensitivity for TQ excitations, are similar to the one found for the 2p 3 3s 5 S o term in O I [13].

4. On a larger Multireference set for 2p 2 ( 3 P )3s 4 P We investigate the reliability of the theoretical parameters by extending the multireference set. Amongst the three terms considered, we focus on the even parity 2p 2 ( 3 P )3s 4 P one, realizing that it is one for which the hyperfine interaction parameters

Comparison with experiments

The hyperfine constants A J and B J are estimated for both isotopes 14 N and 15 N from the hyperfine structure parameters calculated with the most elaborate correlation models (last lines of Tables 2, 3 and 4), using the nuclear data taken from Stone [17] and summarized in Table 5. From equations ( 6)-( 8), one realizes that the ratio between the magnetic hyperfine constants characterizing a given J-level of the two isotopes should -0.28318884 (5) 0.0 be

A J ( 15 N)/A J ( 14 N) = µ I ( 15 N)I( 14 N) µ I ( 14 N)I( 15 N) = -1.4028 . ( 15
)
These non relativistic (NR) theoretical A J and B J hyperfine constants are reported in Table 6 and compared with those derived by Jennerich et al. [1] using saturated absorption spectroscopy ‡. The comparison reveals huge discrepancies between theory and observation. Inconsistencies appear not only in the magnitude of the parameters, but even in the relative signs of the parameters for the different J-levels arising from the same term. ‡ As discussed in [1], the hyperfine constants derived by Cangiano et al. [2] for 4 P J and 4 P o J agree qualitatively with those of Jennerich et al. [1].

Rewriting equations (1)-(4) as

A orb J ≡ G µ µ I I a l K orb J , (16)
A sd J ≡ 1 2 G µ g s µ I I a sd K sd J , (17)
A c J ≡ 1 6 G µ g s µ I I a c K c J , (18)
B J ≡ -G q Q b q K J , (19)
the numerical factors K orb J , K sd J , K c J and K J are calculated from the expectation values of the angular momentum scalar products and reported in Table 7 for each J-values of the (odd and even) 4 P and 4 D terms. The consistency of the theoretical values Table 7. J-dependence of the hyperfine contributions to A J and B J constants for 4 P and 4 D terms. 1 calculated for the different J-values within a given LS term with the relative weights of the different contributions making the total hyperfine constant, is obviously satisfied, by construction. However there is no such constraint in the experimental analysis and one can show that no physical set of underlying {a l , a sd , a c } parameters fits the experimental A J -values with the numerical factors of Table 7.

J K orb J K sd J K c J K J K orb J K sd J K c J K J 4 P 4 D 1/2 -2 3
Some ambiguity was pointed out in the line assignment of the spectra of the 4 P 1/2 → 4 D o 1/2 transition for both isotopes. Due to this identification problem, Jennerich et al [1] determined two possible values for each of the hyperfine constants of the 4 P 1/2 and 4 D o 1/2 states, both reported in Table 6. On the basis of crossover intensity arguments, they gave their preference to the first set. Looking to our theoretical values, we claim that their choice was definitely not the good one. Another surprise appears: they deduced a large B value for 4 D o 3/2 14 N that contrasts with the zero numerical factor K 3/2 of Table 7.

Relativistic corrections

Relativistic effects influence atomic wave functions in basically two ways: through contraction of radial orbitals and through LS term mixing. Contraction effects may be large for high Z, but remain relatively unimportant in light elements like Nitrogen [18]. However, term mixing may affect the wave function and the computed hyperfine interaction constants in significant ways, especially when fine-structure levels belonging to different terms are closely spaced. To investigate term mixing we first perform reference MCHF calculations for configuration expansions generated by SD-excitations from, respectively, {2s 2 2p 2 3s, 2p 4 3s} and {2s 2 2p 2 3p, 2p 4 3p} to active sets (4s3p2d1f ). The resulting non-relativistic radial orbitals are converted to Dirac spinors using the Pauli approximation [19] P (nκ; r) = P MCHF (nl; r) ,

Q(nκ; r) = α 2 d dr + κ r P (nκ; r) . (20)
Here α is the fine structure constant and κ the relativistic quantum number

κ = (j + 1/2)η when l = j + η/2, η = ±1 . (21)
This is followed by relativistic configuration interaction (RCI) calculations for configuration expansions generated by SD-excitations from the multireference sets above to the Dirac spinors (see ref. [20,21] for details about the conversion of radial orbitals and the relativistic computer codes). In relativistic theory only J is a good quantum number and the relativistic configuration expansions account for LS term mixing. The hyperfine interaction constants obtained from the MCHF and RCI wave functions are displayed in Table 8. One observes that the relativistic effects are far from negligible, changing the hyperfine parameters in many cases by more than 10 %. Assuming that the differences between the relativistic and non-relativistic values of the hyperfine constants in the limited calculations with active sets (4s3p2d1f ) are representative of the true differences, we add the former to the non-relativistic results in Table 6 to obtain the final relativistically corrected values (in column "+RC") of the hyperfine interaction constants. It is clear that relativistic effects cannot explain the huge global theoryobservation conflict.

Conclusion

The strong disagreement between theory and observation [1] is really disconcerting. Hyperfine parameters have indeed been estimated ab initio successfully using similar methods for different atomic systems such as Carbon, Nitrogen, Oxygen or Sodium, as described in the introduction. The present calculated hyperfine structure constants disagree so strongly with experiments in comparison of the achieved theoretical convergence of the hyperfine parameters that we presently cast doubt on Cangiano et al's [2] and Jennerich et al's [1] analysis. We therefore encourage further experimental spectroscopic studies and/or reinterpretation of the near-infrared spectra. Last but not least, the extracted isotope shifts from the same spectra revealed a significant unexpected J-dependence of the specific mass shifts in both multiplets. However, the authors themselves [1] observed that the experimental isotope shift values are critically dependent on the correct interpretation of the hyperfine structures of the 14 N and 15 N spectra. The present questioning on the experimental determination of the hyperfine parameters could also be relevant in their isotope shift discussion.

Table 1 .

 1 SD-MR-MCHF and CI total energies, specific mass shift and hyperfine interaction parameters (all in atomic units) as a function of the orbital active set, for 2p 2 3s 4 P . NCSF is the number of CSFs in the wavefunction expansion.

	AS	NCSF	E	S sms	a l	a sd	a c	b q
				SD-MR-MCHF
	HF	1 -54.032 303 -2.5647 3.7993 0.7599 3.6055 1.5197
	3s2p1d	149 -54.116 873 -2.3529 3.5643 0.7310 8.4493 1.2085
	4s3p3d1f	652 -54.165 015 -1.3937 3.7020 0.7999 9.2913 1.2704
	5s4p3d2f1g	1626 -54.189 786 -1.4252 3.6990 0.7837 7.3833 1.2072
	6s5p4d3f2g	3082 -54.196 747 -1.4254 3.6862 0.7713 8.7173 1.2108
	7s6p5d4f3g	5020 -54.199 972 -1.4238 3.6869 0.7723 8.1898 1.2281
	8s7p6d5f3g	7113 -54.201 188 -1.4238 3.6879 0.7751 8.0753 1.2248
	9s8p7d6f3g	9572 -54.201 746 -1.4237 3.6878 0.7752 8.0941 1.2208
	10s9p8d6f3g	11728 -54.201 953 -1.4234 3.6884 0.7743 8.1090 1.2231
				SDTQ-CI
	∪ 4s3p2d1f	40685 -54.203 791 -1.4282 3.6412 0.7668 9.7166 1.1502
	∪ 5s4p3d2f	106472 -54.204 677 -1.4427 3.6170 0.7624 10.511 1.1156
	∪ 6s5p4d3f	210533 -54.204 967 -1.4470 3.6090 0.7609 10.780 1.1049

Table 2 .

 2 SD-MR-MCHF and CI total energies, specific mass shift and hyperfine interaction parameters (all in atomic units) as a function of the orbital active set, for 2p 2 3p 4 P o . NCSF is the number of CSFs in the wavefunction expansion.

	AS	NCSF	E	S sms	a l	a sd	a c	b q
				SD-MR-MCHF
	HF	1 -53.984 055 -2.5786 1.9609 -0.3749 0.0000 -0.7498
	3s2p1d	583 -54.062 052 -2.2945 1.9466 -0.3787 5.1320 -0.7123
	4s3p2d1f	3879 -54.112 044 -1.2732 1.9752 -0.4086 0.5548 -0.6912
	5s4p3d2f1g	10078 -54.134 423 -1.3540 1.9653 -0.3884 0.3990 -0.6708
	6s5p4d3f2g	19200 -54.143 169 -1.3318 1.9673 -0.3839 1.1424 -0.6841
	7s6p5d4f3g	31245 -54.147 455 -1.3170 1.9709 -0.3835 0.7422 -0.6879
	8s7p6d5f3g	44096 -54.148 983 -1.3187 1.9705 -0.3837 0.6826 -0.6887
	9s8p7d6f3g	59110 -54.149 672 -1.3185 1.9702 -0.3850 0.7044 -0.6864
	10s9p8d6f3g	72070 -54.149 921 -1.3188 1.9704 -0.3843 0.7449 -0.6849
				SDTQ-CI
	∪ 4s3p2d1f	124029 -54.150 179 -1.3136 1.9710 -0.3844 0.8878 -0.6848
	∪ 5s4p3d2f	252690 -54.150 330 -1.3139 1.9714 -0.3844 0.9168 -0.6845
	∪ 6s5p4d3f	459494 -54.150 406 -1.3138 1.9717 -0.3843 0.9429 -0.6845
	A J are the most sensitive to TQ excitations due to the large contribution of the contact
	term. The following multireference (MR')
	{2s 2 2p 2 3s, 2p 4 3s, 2s2p 4 , 2s2p 3 3p, 2s2p 2 3s3d}			(14)
	is selected, after a detailed analysis of the eigenvector weights obtained with the first

approach. The results are reported in Table 4. The comparison of the last lines of Table 1 (MR) and Table 4 (MR') illustrates the global stability of the hyperfine parameters, the largest variation (3.6%) being observed for the a c contact contribution.

Table 3 .

 3 SD-MR-MCHF and CI total energies, specific mass shift and hyperfine interaction parameters (all in atomic units) as a function of the orbital active set, for 2p 2 3p 4 D o . NCSF is the number of CSFs in the wavefunction expansion.

	AS	NCSF	E	S sms	a l	a sd	a c	b q
				SD-MR-MCHF
	HF	1 -53.987 051 -2.5789 3.9263 0.7462 0.0000 1.4926
	3s2p1d	720 -54.064 819 -2.2937 3.8979 0.7539 5.1069 1.4194
	4s3p2d1f	4747 -54.114 835 -1.2729 3.9574 0.8137 0.5229 1.3775
	5s4p3d2f1g	12197 -54.137 237 -1.3541 3.9575 0.7769 0.3613 1.3421
	6s5p4d3f2g	23072 -54.146 073 -1.3259 3.9674 0.7698 1.1034 1.3472
	7s6p5d4f3g	37372 -54.150 182 -1.3163 3.9688 0.7661 0.6698 1.3768
	8s7p6d5f3g	52558 -54.151 947 -1.3186 3.9688 0.7658 0.6258 1.3800
	9s8p7d6f3g	70266 -54.152 646 -1.3188 3.9678 0.7692 0.6472 1.3730
	10s9p8d6f3g	85436 -54.152 903 -1.3191 3.9682 0.7688 0.6864 1.3686
				SDTQ-CI
	∪ 4s3p2d1f	156938 -54.153 165 -1.3139 3.9695 0.7691 0.8307 1.3684
	∪ 5s4p3d2f	334710 -54.153 318 -1.3142 3.9706 0.7690 0.8600 1.3678
	∪ 6s5p4d3f	620614 -54.153 406 -1.3138 3.9711 0.7689 0.8848 1.3677

Table 4 .

 4 SD-MR'-MCHF and CI total energies, specific mass shift and hyperfine interaction parameters (all in atomic units) as a function of the orbital active set, for 2p 2 3s 4 P . NCSF is the number of CSFs in the wavefunction expansion.

	AS	NCSF	E	S sms	a l	a sd	a c	b q
				SD-MR'-MCHF
	HF	1 -54.032 303 -2.5647 3.7993 0.7599 3.6055 1.5197
	3s2p1d	461 -54.117 540 -2.3595 3.5468 0.7276 9.0794 1.1821
	4s3p2d1f	2543 -54.168 450 -1.3886 3.6009 0.7798 5.8010 1.1433
	5s4p3d2f1g	6447 -54.192 356 -1.4432 3.6402 0.7761 9.6255 1.1132
	6s5p4d3f2g	12210 -54.199 795 -1.4511 3.6101 0.7584 11.3959 1.0968
	7s6p5d4f3g	19832 -54.203 174 -1.4498 3.6078 0.7589 10.8596 1.1146
	8s7p6d5f3g	28014 -54.204 436 -1.4505 3.6080 0.7615 10.7556 1.1106
	9s8p7d6f3g	37595 -54.205 011 -1.4509 3.6075 0.7615 10.7965 1.1068
	10s9p8d6f3g	45958 -54.205 221 -1.4505 3.6084 0.7607 10.7996 1.1094
				SDTQ-CI
	∪ 4s3p2d1f	93412 -54.205 400 -1.4538 3.5997 0.7592 11.0883 1.0953
	∪ 5s4p3d2f	208694 -54.205 456 -1.4540 3.5977 0.7588 11.1538 1.0922
	∪ 6s5p4d3f	393284 -54.205 469 -1.4543 3.5972 0.7588 11.1705 1.0913

Table 5 .

 5 Nuclear data for 14 N and 15 N[17].

	isotope I	µ I (nm)	Q (b)
	14 N	1 +0.40376100(6) +0.02001(10)
	15 N	1 2

Table 6 .

 6 Hyperfine constant comparison between observation[1] and theory (present work) : non relativistic (NR) ab initio and relativistically corrected (+RC) results. All values are in MHz.

				15 N					14 N
	Levels		exp. [1]	theory		exp. [1]			theory
				NR	+RC				NR		+RC
	4 P 1/2		A +103.4(14) -153.1(23) a	A -139.85	A -140.56		A -69.76(90) +112.3(13) a	B 0.0	A 99.70	B 0.0	A 100.21	B 0.0
	4 P 3/2		-47.93(48)	-88.29	-87.62		35.52(44)	-0.98(48)	62.94	4.10	62.46	4.10
	4 P 5/2		-90.71(71)	-174.75	-175.12		64.76(42)	-3.9(10)	124.58	-5.13	124.84	-5.12
	4 P o 1/2		167.1(13)	75.24	73.29		-133.2(22)	0.0	-53.64	0.0	-52.25	0.0
	4 P o 3/2		70.0(12)	-68.15	-71.60		-48.56(74)	8.69(87)	48.58	-2.58	51.04	-2.95
	4 P o 5/2		46.20(74)	-41.11	-46.52		-32.83(44)	5.0(11)	29.30	3.22	33.16	2.57
	4 D o 1/2		+153.1(23) -103.4(14) a	-106.89	-104.02		-112.3(13) +69.76(90) a	0.0	76.20	0.0	74.15	0.0
	4 D o 3/2		92.4(17)	-49.14	-44.49		-64.41(79)	10.46(88)	35.03	0.0	31.71	0.30
	4 D o 5/2		41.5(14)	-56.74	-51.57		-28.19(62)	-0.2(15)	40.45	-2.30	36.76	-1.69
	4 D o 7/2		-9.35(55)	-77.76	-78.04		6.31(72)	-12.6(13)	55.43	-6.43	55.63	-6.44

a Second proposition of Jennerich et al (see text)

Acknowledgments

TC has a PhD fellowship from the "Fonds pour la formation à la Recherche dans l'Industrie et dans l'Agriculture" of Belgium (Boursier F.R.S. -FNRS). PJ acknowledges support from the Swedish research council. The Fonds de la Recherche Scientifique de Belgique FRFC Convention) and the Communauté française de Belgique (Actions de Recherche Concertées) are acknowledged for their financial support.


Adele Naddeo
email: naddeo@sa.infn.it

Roberta Citro
email: citro@sa.infn.it

Quantum Bose Josephson Junction with binary mixtures of BECs

Keywords: numbers: 03.75.Lm, 67.85.Fg, 74.50.+r

We study the quantum behaviour of a binary mixture of Bose-Einstein condensates (BEC) in a double-well potential starting from a two-mode Bose-Hubbard Hamiltonian. We focus on the small tunneling amplitude regime and apply perturbation theory up to second order. Analytical expressions for the energy eigenvalues and eigenstates are obtained. Then the quantum evolution of the number difference of bosons between the two potential wells is fully investigated for two different initial conditions: completely localized states and coherent spin states. In the first case both the short and the long time dynamics is studied and a rich behaviour is found, ranging from small amplitude oscillations and collapses and revivals to coherent tunneling. In the second case the short-time scale evolution of number difference is determined and a more irregular dynamics is evidenced. Finally, the formation of Schroedinger cat states is considered and shown to affect the momentum distribution.

I. INTRODUCTION

The experimental discovery of Bose-Einstein condensation [1] in dilute systems of trapped alkali-metal atoms, such as rubidium (Rb), lithium (Li), sodium (Na) and ytterbium (Y b), has spurred a renewed interest into the investigation of macroscopic quantum phenomena and interference effects, allowing for a deeper understanding of the conceptual foundations of quantum mechanics [2]. This fascinating research area has been growing up thanks to the high degree of experimental manipulation and control [3]. Interference between condensates released in a potential with a barrier was first observed in 1997 [4] and that paved the way for further investigations on the problem of Bose condensates in a double well potential. Then Josephson oscillations have been observed in one dimensional optical potential arrays [5]. A single bosonic Josephson junction was produced for the first time in 2005 with Rb atoms and its dynamics was experimentally investigated both within tunneling as well as self-trapping regime [6][7] [8]. More recently, mixtures of 85 Rb and 87 Rb atoms have been produced and experimentally investigated [9] as well, whose intraspecies scattering lengths could be tunable via magnetic and optical Feshbach resonances. Furthermore the realization of heteronuclear mixtures of 87 Rb and 41 K atoms with tunable interspecies interactions [10] paved the way to the exploration of double species Mott insulators and, in general, of the quantum phase diagram of two species Bose-Hubbard model [11]. The interplay between the interspecies and intraspecies scattering produces deep consequences on the properties of the condensates, such as the density profile [12] and the collective excitations [13]. However, the wide tunability of such interactions makes a BEC mixture a very interesting subject of investigation, both from experimental and theoretical side as a mean of studying new macroscopic quantum tunneling phenomena as well as the interplay between quantum coherence and nonlinearity. Indeed novel and richer behaviours are expected in such a multicomponent BEC.

On the theoretical side, a bosonic Josephson junction with a single species of BEC has been widely investigated by means of a two-mode approximation [14][15] [16], within the classical as well as the quantum regime. In the classical regime, characterized by large particle numbers and weak repulsive interactions, the Gross-Pitaevskii equation provides a reliable description. Within the two mode approximation it reduces to two generalized Josephson equations which describe the time evolution of the relative phase and the population imbalance between the wells [15] and differ from their superconducting counterpart [START_REF] Barone | Physics and Applications of the Josephson Effect[END_REF] by the presence of a nonlinear term which couples the variables. Because of such a term, a bosonic Josephson junction exhibit a variety of novel phenomena which range from π-oscillations to macroscopic quantum self-trapping (MQST) [15]. While the π-oscillations, as well the usual Josephosn ones, deal with a symmetric oscillation of the condensate about the two wells, the MQST phenomenon is characterized by a broken symmetry phase with a population imbalance between the wells. In the quantum regime, characterized by smaller values of the particle number and strong interactions, an increasing of phase fluctuations is observed together with the suppression of number fluctuations. Furthermore the time evolution is characterized by phase collapse and revival [START_REF] Javanainen | [END_REF]. The quantum behaviour of bosonic Josephson junctions has been deeply investigated by means of the usual quantum phase model [19] [20] [21] as well as by starting from a two-mode Bose-Hubbard Hamiltonian [22][23] [24]. In this context the phase coherence of the junction has been characterized by studying the momentum distribution [25] [20]. The generation and detection of Schroedinger cat states has been investigated as well; indeed the presence of such a kind of states reflects in the strong reduction of the momentum-distribution contrast [24] [26].

More recently such a theoretical analysis has been successfully extended to a binary mixture of BECs in a double well potential [27][28] [29] [30]. The semiclassical regime in which the fluctuations around the mean values are small has been deeply investigated and found to be described by two coupled Gross-Pitaevskii equations. By means of a two-mode approximation such equations can be cast in the form of four coupled nonlinear ordinary differential equations for the population imbalance and the relative phase of each species. The solution results in a richer tunneling dynamics. In particular, two different MQST states with broken symmetry have been found [29], where the two species localize in the two different wells giving rise to a phase separation or coexist in the same well respectively. Indeed, upon a variation of some parameters or initial conditions, the phase-separated MQST states evolve towards a symmetry-restoring phase where the two components swap places between the two wells, so avoiding each other. Recently, the coherent dynamics of a two species BEC in a double well has been analyzed as well focussing on the case where the two species are two hyperfine states of the same alkali metal [31].

In this paper we study the quantum behaviour of a binary mixture of Bose-Einstein condensates (BEC) in a double-well potential starting from a two-mode Bose-Hubbard Hamiltonian. We analyze in detail the small tunneling amplitude regime where number fluctuations are suppressed and a Mott-insulator behaviour is established. We perform a perturbative calculation up to second order in the tunneling amplitude and study the stationary states and the dynamics of the two species bosonic Josephson junction. Finally, the dynamical generation of Schroedinger cat states is investigated starting from an initial coherent spin state and shown to affect the time-dependent population imbalance and momentum distribution [25] [20]. We focus on the contrast in the momentum distribution between the two wells and show how it vanishes for a two-component cat state. That could be interesting in view of the experimental realization of macroscopic superpositions of quantum states [26] [32].

The paper is organized as follows. In Section 2 we introduce our model Hamiltonian within the two-mode approximation and define the various parameters. Then we adopt the angular momentum representation and focus on the small tunneling amplitude regime. In Section 3 we apply perturbation theory in the tunneling amplitude to our Hamiltonian and find analytical expressions for the energy eigenvalues and eigenstates up to second order. Section 4 and 5 are devoted to the study of the quantum evolution of the number difference of bosons between the two wells in correspondence of two different initial conditions: completely localized states and coherent spin states. In the first case both the short and the long time dynamics is studied and a rich behaviour is evidenced, ranging from small amplitude oscillations and collapses and revivals to coherent tunneling. In the second case the short-time scale evolution of number difference is determined and a more irregular dynamics is evidenced, with suppression of the dominant frequency when the number of bosons increase. Then, Schroedinger cat states are shown to generate as a result of the time-evolution of an initial coherent state when the tunneling between the two wells is suppressed, and their influence on the contrast in the momentum distribution is studied. Finally, in Section 6 some conclusions and outlooks of this work are presented.

II. THE MODEL

A binary mixture of Bose-Einstein condensates [28][29] loaded in a double-well potential is described by the general many-body Hamiltonian:

H = H a + H b + H ab (1)
where

H a = d -→ r - h2 2m a ψ + a ∇ 2 ψ a + ψ + a V a ( -→ r ) ψ a + 1 2 d -→ r d - → r ψ + a ( - → r ) ψ + a -→ r U aa -→ r -- → r ψ a - → r ψ a ( - → r ) , (2)
H b = d -→ r - h2 2m b ψ + b ∇ 2 ψ b + ψ + b V b ( - → r ) ψ b + 1 2 d - → r d - → r ψ + b ( - → r ) ψ + b -→ r U bb - → r -- → r ψ b - → r ψ b ( - → r ) (3)
are the Hamiltonians for bosons of species a and b respectively and

H ab = d - → r d -→ r ψ + a ( - → r ) ψ + b -→ r U ab - → r --→ r ψ a -→ r ψ b ( - → r ) , (4)
is the interaction term between bosons of different species. For dilute mixtures one can replace the interaction potentials U aa , U bb and U ab with the effective contact interactions:

U aa -→ r --→ r = g aa δ - → r --→ r , U bb - → r -- → r = g bb δ -→ r --→ r , U ab - → r -- → r = g ab δ - → r -- → r , (5)
where g aa = 4πh 2 aaa ma and g bb = 4πh 2 a bb m b are the intraspecies coupling constants of the species a and b respectively, m a and m b being the atomic masses and a aa , a bb the s-wave scattering lengths; furthermore g ab = 2πh 2 a ab m ab is the interspecies coupling constant, where m ab = mam b ma+m b is the reduced mass and a ab is the associated s-wave scattering length. In this way the Hamiltonian (1)-( 4) can be rewritten as:

H i = d - → r - h2 2m i ψ + i ∇ 2 ψ i + ψ + i V i ( -→ r ) ψ i + g ii 2 d - → r ψ + i ψ + i ψ i ψ i ; i = a, b (6)
H ab = g ab d - → r ψ + a ψ + b ψ a ψ b . (7)
Here V i ( -→ r ) is the double well trapping potential and, in the following, we assume

V a ( - → r ) = V b ( -→ r ) = V ( - → r ); ψ + i ( - → r ) , ψ i ( -→ r ), i = a,
b are the bosonic creation and annihilation operators for the two species, which satisfy the commutation rules:

ψ i ( -→ r ) , ψ j -→ r = ψ + i ( - → r ) , ψ + j - → r = 0, (8)
ψ i ( -→ r ) , ψ + j -→ r = δ ij δ -→ r -- → r , i,j = a, b, (9)
and the normalization conditions:

d -→ r |ψ i ( -→ r )| 2 = N i ; i = a, b, (10)
N i , i = a, b being the number of atoms of species a and b respectively. The total number of atoms of the mixture is N = N a + N b . Furthermore in the following we refer to values of the coupling constants g aa , g bb and g ab such that the condition of stability of the binary mixture, g aa g bb ≥ g 2 ab , is met and the phase separation regime is avoided. Now a weak link between the two wells produces a small energy splitting between the mean-field ground state and the first excited state of the double well potential and that allows us to reduce the dimension of the Hilbert space of the initial many-body problem. Indeed for low energy excitations and low temperatures it is possible to consider only such two states and neglect the contribution from the higher ones, the so called two-mode approximation [14] [15] [16]. In this way, by taking into account for each of the two species a and b the mean-field ground states φ a g , φ b g and the mean-field excited states φ a e , φ b e , the wave functions ψ i , i = a, b, can be rewritten as:

ψ a = a g φ a g + a e φ a e ψ b = b g φ b g + b e φ b e , (11)
where d -→ r φ i g(e)

2

= 1, i = a, b, and a + g , b + g and a + e , b + e (a g , b g and a e , b e ) are the creation (annihilation) operators for a particle of the species a, b in the ground and the excited state respectively. They satisfy the usual bosonic commutation relations a i , a

+ j = b i , b + j = δ ij . Furthermore φ i
g,e , i = a, b are assumed real for simplicity and such that φ i g 3 φ j e = φ i e 3 φ j g = 0, which simplifies the calculations. Let us change the basis and switch to the atom number states in such a way that the expectation value of the population of the left and right well can be defined.

The new annihilation operators are a

L = 1 √ 2 (a g + a e ), a R = 1 √ 2 (a g -a e ) and b L = 1 √ 2 (b g + b e ), b R = 1 √ 2 (b g -b e )
for the species a and b respectively, so that the wave functions (11) become:

ψ a = 1 √ 2 a L φ a g + φ a e + 1 √ 2 a R φ a g -φ a e ψ b = 1 √ 2 b L φ b g + φ b e + 1 √ 2 b R φ b g -φ b e . ( 12
)
By substituting Equations ( 12) into the Hamiltonian ( 6)- (7), after some algebra we obtain its second quantized version within the two-mode approximation:

H = E a c 8 a + R a R -a + L a L 2 - E a J N a a + R a L + a + L a R + δE a a + R a L + a + L a R 2 + E b c 8 b + R b R -b + L b L - E b J N b b + R b L + b + L b R + δE b b + R b L + b + L b R 2 + 1 4 Λ ab a + L a L -a + R a R b + L b L -b + R b R + - 1 4 a + R a L + a + L a R b + R b L + b + L b R κ ab e,g + κ ab g,e -κ ab g,g -κ ab e,e + 1 2 N a E a g + E a e + ( 1 3 ) + 1 4 N a (N a -2) κ a g,g + κ a e,e + N 2 aL + N 2 aR -N a κ a g,e + 1 2 N b E b g + E b e + + 1 4 N b (N b -2) κ b g,g + κ b e,e + N 2 bL + N 2 bR -N b κ b g,e + 1 4 N a N b κ ab e,g + κ ab g,e + κ ab g,g + κ ab e,e ,
where

N i = N iL + N iR , i = a, b
, is the number of atoms of species a and b respectively, expressed as a sum of numbers of atoms in the left and right well. The parameters are defined as follows:

E i g = d - → r - h2 2m i φ i g ∇ 2 φ i g + φ i g V φ i g ; i = a, b (14)
E i e = d - → r - h2 2m i φ i e ∇ 2 φ i e + φ i e V φ i e ; i = a, b (15)
κ a i,j = g aa 2 d - → r |φ a i | 2 φ a j 2 ; i, j = g, e (16)
κ b i,j = g bb 2 d - → r φ b i 2 φ b j 2 ; i, j = g, e (17)
κ ab i,j = g ab d -→ r |φ a i | 2 φ b j 2 ; i, j = g, e (18)
Λ ab = 4g ab d - → r φ a g φ a e φ b g φ b e ; ( 1 9 )
E i c = 4κ i g,e ; i = a, b (20)
δE i = κ i g,g + κ i e,e -2κ i g,e 4 ; i = a, b (21)
E i J = N i 2 E i e -E i g + N i 2 (N i -1) κ i g,g -κ i e,e + N j 2 κ ab e,e -κ ab g,g + κ ab g,e -κ ab e,g ; i, j = a, b; i = j. (22)
In Eq. ( 13), the terms proportional to E i J , i = a, b, describe tunneling of particles of species a and b from one to the other well while the terms proportional to E i c , i = a, b, deal with the local interaction within the two wells and the terms proportional to δE i correspond to additional two-particle processes. Finally the terms proportional to Λ ab and κ ab i,j couple the two species and then various constant terms follow, which we will drop for simplicity.

In this paper we focus on the small tunneling amplitude regime where number fluctuations are suppressed and a Mott-insulator behaviour is established, so it is convenient to introduce the angular momentum representation for the species a and b as follows:

J a x = 1 2 a + R a L + a + L a R , J a y = i 2 a + R a L -a + L a R , J a z = 1 2 a + R a R -a + L a L , J b x = 1 2 b + R b L + b + L b R , J b y = i 2 b + R b L -b + L b R J b z = 1 2 b + R b R -b + L b L , (23)
where the operators J a i , J b i , i = x, y, z, obey to the usual angular momentum algebra and the following relations hold:

(J a ) 2 = Na 2 Na 2 + 1 , J b 2 = N b 2 N b 2 + 1 . (24)
In particular, the components

J i z = 1 2 (N iR -N iL ), i = a, b
give the difference of the number of bosons of the species i, N iL and N iR , occupying the two minima of the double well potential, i. e. the population imbalances, which are experimentally observable quantities. Thus Hamiltonian ( 13) can be cast in the following form:

H = E a c 2 (J a z ) 2 -2 E a J N a J a x + 4δE a (J a x ) 2 + E b c 2 J b z 2 -2 E b J N b J b x + 4δE b J b x 2 + +Λ ab J a z J b z -J a x J b x κ ab e,g + κ ab g,e -κ ab g,g -κ ab e,e , (25)
where the constant terms have been dropped. Let us now simplify the notation by introducing the following parameters:

Λ a = E a c , C a = 4δE a , K a = 2 E a J Na , Λ b = E b c , C b = 4δE b , K b = 2 E b J N b , D ab = κ ab e,g + κ ab g,e -κ ab g,g -κ ab e,e (26)
and rewrite the Hamiltonian (25) as:

H = 1 2 Λ a (J a z ) 2 -K a J a x + C a (J a x ) 2 + 1 2 Λ b J b z 2 -K b J b x + C b J b x 2 + +Λ ab J a z J b z -D ab J a x J b x . (27)
Within the experimental parameters range it is possible to show that C i Λ i , K i , i = a, b, and D ab Λ ab [7][29], then in the following we put C a = C b = 0 and D ab = 0, which corresponds to neglecting the spatial overlap integrals between the localized modes in the two wells. In this way the binary mixture of BECs within a simple two-mode approximation maps to two Ising-type spin model in a transverse magnetic field. It is possible, also in this case of binary mixtures of BECs, to take into account the overlap between wave functions localized in different wells and then retain the terms just neglected, as done in Ref. [16] for the single species Bose Josephson junction. In this way an improved two-mode model can be defined, which gives rise to a quantitative agreement with the results of a numerical integration of the coupled Gross-Pitaevskii equations describing the system [33]. We would point out that throughout the paper we adopt a simple two-mode approximation because we are interested in the qualitative features of the quantum dynamics of a binary mixture of BECs in a double well.

In the following we will focus on the symmetric case Λ a = Λ b = Λ and K a = K b = K because it allows us to perform analytical calculations while capturing many relevant phenomena characterizing the physics of the system. An investigation of the full parameter space will require numerical calculations and will be the subject of a future publication together with a detailed comparison between analytical and numerical results and a quantitative evaluation of the validity of the two-mode approximation [34]. So the model Hamiltonian (25) becomes:

H = H 0 + H I , (28)
H 0 = 1 2 Λ (J a z ) 2 + 1 2 Λ J b z 2 + Λ ab J a z J b z , (29)
H I = -K J a x + J b x , (30)
where, in the small tunneling amplitude regime, H I is considered as a perturbation. The total Hamiltonian commutes with (J a ) 2 and J b 2 , which leads to the conservation of total angular momentum with quantum numbers j a = Na

b = -N b 2 , ..., N b 2 .
As a first step we need to diagonalize the unperturbed Hamiltonian (29), which can be done by performing the following θ = π 4 rotation on the operators J a z , J b z :

O 1 z = a 1 J a z -a 2 J b z J 2 z = a 1 J a z + a 2 J b z , a 1 = a 2 = 1 √ 2 , (31)
while an analogous rotation needs to be carried out on J a x , J b x entering the perturbation (30). As a result we get:

H = 1 2 (Λ -Λ ab ) O 1 z 2 + 1 2 (Λ + Λ ab ) J 2 z 2 - 2K √ 2 J 2 x , (32)
which, by defining Λ 1 = Λ -Λ ab , Λ 2 = Λ + Λ ab , and

O 1 z = O 1 z √ 2 , O 1 x = O 1 x √ 2 , J 2 z = J 2 z √ 2 , J 2 x = J 2 x √ 2
, can be cast in the final form:

H = Λ 1 O 1 z 2 + Λ 2 J 2 z 2 -2K J 2 x . (33)
In the following Section we will find analytical expressions for the eigenvalues and the eigenvectors up to second order by performing perturbation theory in the tunneling amplitude.

III. STATIONARY STATES

In the present Section we apply second-order perturbation theory to the Hamiltonian of Eq. ( 33) in the small tunneling amplitude limit, which allows us to derive analytical expressions for the stationary states of the system.

In order to pursue this task let us rewrite Eq. ( 33) in dimensionless form by assuming Λ1 2 as unit of energy:

H = 2 O 1 z 2 + 2λ J 2 z 2 -2k J 2 x , (34)
where λ = Λ2 Λ1 and k = 2K Λ1 , then take

H 0 = 2 O 1 z 2 + 2λ J 2 z 2 (35)
as unperturbed Hamiltonian and

H I = -2k J 2 x , (36)
as a small perturbation term. Here J 2 i , i = x, y, z, obey the usual angular momentum algebra and the following relation holds:

J 2 2 = N 2 2 N 2 2 + 1 , (37)
where:

N 2 = N a + N b 2 . ( 38
)
In principle, the rotated basis

{|m 1 , m 2 } = m 1 = 1 2 (m a -m b ) , m 2 = 1 2 (m a + m b ) of the unperturbed Hamil- tonian (35) is given by the eigenvectors of O 1 z = J a z -J b z 2 ( O 1 z |m 1 = m 1 |m 1 ) and J 2 z ( J 2 z |m 2 = m 2 |m 2 ) with m 1 = -|ja-j b | 2 , ..., |ja-j b | 2 and m 2 = -(ja+j b ) 2 , ..., (ja+j b ) 2 , whose corresponding eigenvalues are E (0) m1,m2 = 2 (m 1 ) 2 + 2λ (m 2 ) 2 .
The presence of the operator O 1 z , which does not commute with the perturbation term H I , makes the problem of finding eigenvalues and eigenvectors of the full Hamiltonian (34) within perturbation theory much more involved. In order to simplify the treatment and carry out analytical calculations while retaining the relevant phenomenology, we concentrate on the particular case of a binary mixture where the two species are equally populated, i. e. N a = N b , and have the same population imbalance between the two wells, i. e. m a = m b . This situation allows us to describe the quantum dynamics of the system in correspondence of the MQST regime, for which we need a completely localized initial state. That fixes m 1 = 0 while m 2 = m a = -Na 2 , ..., Na 2 could be an even or odd integer depending on N a even or odd, and leads to the following zero-order eigenvalues:

E (0) 0,m ± 2 = 2λ (m 2 ) 2 = λ 2 (m a + m b ) 2 .
Each eigenvalue is two-fold degenerate, with the only exception of the ground state for N a even, E

0,0 = 0, which is nondegenerate. The two-dimensional subspace of degeneracy is spanned by the states |0, ±m 2 (where J 2 z |0, ±m 2 = ±m 2 |0, ±m 2 ) and the corresponding zero-order eigenvectors are:

h (0) 0,m ± 2 = 0, m ± 2 = 1 √ 2 (|0, m 2 ± |0, -m 2 ) . ( 39
)
By switching on the perturbation term [START_REF] Arecchi | [END_REF] it is possible to show that the degeneration is lifted starting from the levels with smaller m 2 ; in general the double degeneracy of the zero-order eigenvalues E (0) 0,m ± 2 will be lifted at the 2m 2 -th order of perturbation theory [22]. By applying perturbation theory [START_REF] Cohen-Tannoudji | Quantum Mechanics[END_REF] up to order k 2 , we obtain the following corrected eigenvalues:

E (2) 0,m ± 2 = 2λ (m 2 ) 2 + k 2 λ j 2 (j 2 + 1) + (m 2 ) 2 4 (m 2 ) 2 -1 ; m 2 = 1, 1 2 , (40)
E (2) 0,1 ± = 2λ + k 2 λ j 2 (j 2 + 1) + 1 3 ± j 2 (j 2 + 1) 2 ; N 2 even, (41)
E (2) 0, 1 2 ± = λ 2 ∓ k j 2 (j 2 + 1) + 1 4 - k 2 4λ j 2 (j 2 + 1) - 3 4 ; N 2 odd, ( 42
)
where j 2 = N2 2 . Furthermore, for N 2 even, the nondegenerate ground state |0, 0 belongs to the symmetry class of 0, m + 2 . The corresponding eigenvectors, up to order k 2 , are given in the Appendix. As a final remark we would make some qualitative considerations on the precision of the second order perturbative results we just obtained. In the previous Section we pointed out that the tunneling amplitude K (and, then, the renormalized tunneling amplitude k) must be chosen sufficiently small in order to realize a weak-link. Furthermore, we are considering a regime very close to the Mott-insulator one; so, being the onsite energies much larger than the tunneling amplitudes, a perturbation theory in k gives reliable results. In general, the perturbative corrections to the eigenvalues in Eqs. ( 40)-(42) behave as powers of kj 2 ∼ kN 2 for large N 2 . As a consequence the bigger the value of N 2 is, smaller is the value of k at which the perturbative expansion is reliable, as shown in detail in Ref. [22] for the single species Bose-Hubbard dimer. For the ground state, and in general for the lower-energy states, the divergence of the perturbative results abruptly increases with N 2 . So, for BEC mixtures such as the one of JILA setup with N 2 ∼ 6.5 × 10 4 [9] the perturbative results for the lower-energy eigenvalues should be valid for very small values of k. Conversely, for the higher-energy levels (i. e. m 2 = N2

2 ) the second order corrections to Eq. ( 40) go as k 2 N2 N2-1 instead of (kN 2 ) 2 of the low-lying states and then decrease as N 2 increases. As a consequence, the higher order corrections are much more strongly decreasing functions of N 2 and the perturbative results are accurate to larger values of k.

In the following Sections we use the analytical expressions of energy eigenvectors derived in the Appendix, see Eqs. (62)-(69), in order to study the quantum evolution of J 2 z (τ ) = 1 2 J a z (τ ) + J b z (τ ) , that is the number difference of bosons of species a and b between the two wells of the potential.

IV. DYNAMICS: COMPLETELY LOCALIZED INITIAL STATES

In this Section we investigate the quantum evolution of the number difference of bosons of species a and b between the two wells assuming a completely localized state as initial condition. That could be interesting in order to elucidate the quantum behavior of the system in correspondence of the classical MQST regime and to put in evidence new phenomena including quantum coherence in a multicomponent system. In such a case we will study both the short and the long time dynamics: as a result a rich behaviour emerges, ranging from small amplitude oscillations and collapses and revivals to coherent tunneling. Although such a physics is well known for the single component Bose Josephson junction, in our case the dynamics shows that the two species can coexist in the same potential well as if there would be an attractive interaction between them.

As a first step let us recall the general formula which gives the time evolution of the mean value of

J 2 z = 1 2 J a z + J b z [35]: J 2 z (τ ) = n=m ± 2 n =m ± 2 φ * n φ n h 0,n J 2 z h 0,n e i E0,n-E 0,n τ , (43)
where τ = Λ1 2h t is the dimensionless time, the sums are over all the eigenvectors h 0,m ±

2

, being m 2 = 0 or 1 2 , ..., N2 2 , and φ n are the projections of the initial state |ψ (0) on the basis h 0,m ± 2 :

|ψ (0) = n=m ± 2 φ n h 0,n . ( 44
)
So it is clear how the knowledge of eigenvalues and eigenvectors is enough in order to study the quantum evolution of J 2 z , the Bohr frequencies involved, E 0,n -E 0,n , and the corresponding weights φ * n φ n h 0,n J 2 z h 0,n . Let us now study the dynamics of the system when all the bosons of species a and b are initially contained in one of the two wells of the potential, say the right one, and then the imbalances of the two species coincide, so that N aR = N a , N aL = 0, N bR = N b , N bL = 0; furthermore the two species are equally populated, i.e. N a = N b . That implies m 1 = 0 and m 2 = Na 2 = N2 2 in our center of mass rotated basis. The corresponding initial condition is:

|ψ (0) = 0, N 2 2 . ( 45
)
In order to investigate the short timescales evolution we need to keep terms up to second order in the tunneling amplitude k when we compute the weights in Eq. ( 43). We find that:

J 2 z (2) (τ ) = N 2 2 + k 2 N 2 2λ 2 (N 2 -1) 2 [cos(ω µ τ ) -1] , (46)
where the frequency involved is:

ω µ = E (2) 0, N 2 2 ± -E (2) 0,( N 2 2 -1) ± = 2λ (N 2 -1) - k 2 λ N 2 + 1 (N 2 ) 2 -4N 2 + 3 . ( 47
)
At short timescales small amplitude oscillations with frequency ω µ around the initial condition (N 2R = N 2 , N 2L = 0) are observed and that coincides with a strongly self-trapped regime.

In order to investigate the dynamics at longer timescales we have to take into account also the small splittings ∆ E 0,

N 2 2 ± and ∆ E 0,( N 2 2 -1)
± of the two higher pairs of quasidegenerate eigenvalues which provide two further frequencies (see Ref. [23] for the derivation):

ω 0 = ∆ E 0, N 2 2 ± = k N 2 λ N 2 -1 N2 2 N 2 -2 (N2-1)! ω 1 = ∆ E 0,( N 2 2 -1) ± = k N 2 -2 λ N 2 -3 (N2-1)(N2-2) 2 N 2 -4 (N2-3)! . ( 48
)
The whole result is:

J 2 z (2) (τ ) = N 2 2 cos(ω 0 τ ) + k 2 N 2 4λ 2 (N 2 -1) 2 N 2 2 [cos(ω 1 τ ) -cos(ω 0 τ )] +2 cos(ω µ τ ) cos( ω 1 2 τ ) -cos(ω 1 τ ) -cos(ω 0 τ ) , (49)
and, by putting ω 0 = ω 1 = 0, the short timescale dynamics, Eq. ( 46), is recovered. Summarizing, at longer timescales the two species bosons are still localized in the initial potential well but the quantum dynamics exhibits collapses and complete revivals. Indeed the coefficient cos ω1 2 τ , which multiplies the higher frequency term cos (ω µ τ ), gives rise to the beat, which is responsible for the observed collapses and revivals at timescales fixed by ω 1 , as shown in Fig. 1. Finally, at very large timescales determined by the frequency ω 0 all the bosons tunnel coherently back and forth between the two traps; only the first term cos(ω 0 τ ) is responsible of such a coherent tunneling, since all harmonic functions containing the frequency ω 1 and ω µ are small in amplitude and proportional to k 2 , thus they are unable to transfer bosons from one trap to the other.

The tunneling dynamics within macroscopic quantum self-trapping regime described above is analogous to that of the π-mode fixed point obtained by the Gross-Pitaevskii approach [29], where the two species localize in the same well despite the repulsive interaction between them. Let us finally note that, despite the explicit dependence on λ of the frequencies (47)-( 48), the different physics related to the three time scales described above is simply due to the energy splitting introduced by the renormalized tunneling for small Λ ab . Thus in the case of a mixture of BECs with equal population the dynamics remains similar to that of a single component BEC, apart the coexistence of the two species in the same well. This is shown in Fig. 1 where the tunneling dynamics for the single species BEC is also reported.

As for the experimental detection of the long timescales phenomena (collapses/revivals and coherent tunneling), since the time for their appearance is abruptly increased with N 2 , this implies a rapid decrease of the characteristic frequencies rendering more difficult the observation of the intermediate and long time behavior in current BEC experiments. Indeed, by a comparison of Eqs. ( 47)-(48), we clearly see that the frequency ω µ is of order k 0 while ω 0 and ω 1 are of order k N2 and k N2-2 respectively. That reflects on the dimensionless time-scales for such phenomena, which range from order 1 for short time dynamics to 10 7 and 10 11 respectively for the intermediate and long time dynamics. Indeed pure condensates consisting of 1150 ± 150 atoms of 87 Rb loaded in a double well have been recently realized [6][7] thus rendering the detection of the short time behavior possible. Mixtures with a number of atoms ranging from 9 × 10 3 and 5 × 10 3 ( 87 Rb and 41 K [10]) to 4 × 10 4 and 9 × 10 4 ( 85 Rb and 87 Rb [9]) have also been recently realized, but in this case very small characteristic frequencies are implied. However, these phenomena may be relevant for molecular systems where the number of vibrational excited quanta is small, allowing for the experimental observation of intermediate and long time dynamics also.

In the next Section we will further investigate the dynamics of the system by assuming as initial state a simple coherent state and then study the formation of a particular superposition of such coherent states, the so called Schroedinger cat states.

V. DYNAMICS: COHERENT SPIN INITIAL STATES AND SCHROEDINGER CAT STATES

In this Section we choose as initial condition a simple coherent spin state [START_REF] Arecchi | [END_REF] and study the short-time scale evolution of number difference; in this way a more complex dynamics will appear. Finally, we study the generation of Schroedinger cat states; in particular, we focus on the contrast in the momentum distribution and show how it vanishes for a two-component cat state.

Let us start by considering as initial condition the following coherent spin state [START_REF] Arecchi | [END_REF]:

|ψ (0) = C N2/2 m2=-N2/2 N 2 ! N2 2 + m 2 ! N2 2 -m 2 ! tan m2 θ 2 e -im2φ |0, m 2 , (50)
where the coefficient C is:

C = sin N2/2 θ 2 cos N2/2 θ 2 e -i(N2/2)φ , (51)
and θ and φ are two angles characterizing the superposition. The time evolution of the mean value of J 2 z up to first order in the tunneling amplitude k is given by:

J 2 z (1) (τ ) = J 2 z (0) (τ ) + k λ sin(θ) 2 N2 [C 1 [cos(ω e τ ) -1] + C 2 sin (ω e τ ) + N2/2-1 n=0or1/2 N 2 ! N2 2 + n ! N2 2 -n ! N 2 -2n 2 (2n + 1) A n   , (52)
where

A n = tan 2n+1 θ 2 [cos (F n τ + φ) -cos (φ)] - 1 tan 2n+1 θ 2 [cos (F n τ -φ) -cos (φ)] (53)
with frequencies

F n = E (0) 0,(n+1) ± -E (0) 0,(n) ∓ = λ (4n + 2)
. Furthermore the coefficients C 1 and C 2 are given by:

C 1 = N 2 ! N2 2 + 1 ! N2 2 -1 ! cos (φ) N 2 6 - 1 3 tan 3 θ 2 - 1 tan 3 θ 2 - N 2 2 + 1 tan θ 2 - 1 tan θ 2 , (54)
C 2 = N 2 ! N2 2 + 1 ! N2 2 -1 ! tan θ 2 + 1 tan θ 2 N 2 6 - 1 3 sin (3φ) - N 2 2 + 1 sin (φ) , (55)
for N 2 even, and

C 1 = N 2 ! tan 2 θ 2 - 1 tan 2 θ 2 [cos(ω e τ ) -1] + 2 sin (2φ) sin (ω e τ ) , (58)
where the dominant frequency ω e is equal to ω e = E

(2)

0,1 + -E (2) 0,1 -= k 2 λ N2 2 N2 2 + 1 .
The corresponding expression for N 2 odd is: where ω e = E

J 2 z (0) (τ ) = - N 2 2 cos (θ) + 1 2 sin(θ) 2 N2 N 2 ! N2 2 + 1 2 ! N2 2 -1 2 ! tan θ 2 - 1 tan θ 2 [cos(ω e τ ) -1] -2 sin (φ) sin (ω e τ ) , (59)
(2) 0, 1

--E (2) 0, 1 2 + = 2k N2 2 N2 2 + 1 + 1 4 . 2
As one can see, the dominant frequency is gradually suppressed with the number of bosons N 2 = 1 2 (N a + N b ) [22]. This is clearly seen in Fig. 2 where the boson number difference between the two traps is plotted for different N 2 (even) values as a function of the dimensionless time τ . One also notices a decrease of the oscillation amplitude at increasing N 2 . The effect of λ is instead shown in Fig. 3 where the short-time dynamics of the boson number difference is analyzed for two values of the interspecies interaction. When λ increases the amplitude of the oscillations decreases. The detection of the mixture dynamics is thus more favorable for values of λ smaller than unity.

. .

A. Cat states

Let us consider the coherent spin state (50); the expectation value of the Hamiltonian (34) on such state is given by:

ψ (0) | H |ψ (0) = 2λn 2 /2 -2k (N 2 /2) 2 -n 2 cos φ, (60)
where n = -(N 2 /2) cos θ and has the maximum value for φ = 0, θ = π/2. This result also corresponds to the mean-field result for the energy. Now, starting from the coherent spin state (50) we are interested in looking for Schroedinger cat states. Such states are quantum superposition of macroscopic states and their realization has already been suggested for a single species Bose-Josephson junction in [24]. Also in the case of a Bose-Josephson junction with binary mixtures one might realize cat states from the time-evolution of an initially coherent state following a sudden rise of the barrier between the two wells. Thus we consider at time t = 0 + a zero inter-well coupling k, i.e. the time evolution is governed by the Hamiltonian H 0 in Eq. [START_REF] Cohen-Tannoudji | Quantum Mechanics[END_REF]. For each basis vector |0, m 2 of the coherent state (50), the time-evolution is given by |0, m 2 (t) = e -i2πm 2 2 t/T |0, m 2 , where T = hπ/λ is the so-called revival time such that |ψ(T ) = |ψ(0) . Considering now the times T /2p, p integer, the time evolution of the coherent state is governed by the factor exp(-iπm 2 2 /p) which satisfies the property exp(-iπ(m 2 + p) 2 /p) = (-1) p exp(-iπm 2 2 /p), depending on the parity of p. For the choice of even p, a discrete Fourier transform leads to the cat state:

|ψ(T/2p) = p-1 k=0 u k e iπkN2/p |e -i2πk/p ψ , (61)
i.e. a superposition of p coherent states, where u k = 1/p p-1 m2=0 e -iπm 2 2 /p e i2πkm/p . In particular, the cat state affects the momentum distribution. This dependence could be important to probe experimentally their existence. In particular, when considering the two-component cat state, i.e. for the choice p = 2, one obtains that the contrast in the momentum distribution, i.e. the expectation value of J 2

x on the unperturbed state, vanishes [24]. Furthermore, the amplitude of the intervals of time in which the contrast is zero increases with increasing N 2 as clearly shown in Fig. 4.

It should be noted that despite the close similarity in the behavior of the contrast between the single component BJJ and the double one, the mixture could be a better candidate for the creation and detection of cat states. In fact their creation time is πh/λ and since for repulsive interaction between the two species and Λ > Λ ab we get λ = -[(1 + Λ/Λ ab )/(1 -Λ/Λ ab )] > 1, such time can be made short enough to render their detection more favorable. But, in order to establish if the mixture will be an useful candidate for the creation and detection of cat states, it would be necessary to study the influence of the various decoherence mechanisms. In general, the main decoherence mechanism which limits the lifetime and the size of BECs is the three-body recombination. Three-body losses strongly influence the formation and evolution of complex superpositions of states of different phase, i. e. the mechanism of quantum collapse and revival of the phase [37] which is adopted as a protocol to realize cat states [32]. In order to establish the robustness of Schroedinger cat states against three body losses it is then mandatory to compare the decoherence time and the revival time keeping in mind that, if the decoherence time scale is much faster than the revival time, revivals will not occur and intermediate superpositions of states are destroyed. In general, for atomic three-body decay the recombination loss at each site is characterized by the decay rate γ = K3 12 dx |w(x)| 6 , where w(x) is the localized wave function at the site x and K 3 is the three-body rate constant [37]. It can be shown [38] that, for a binary mixture of BECs, such a decay rate takes the form γ BM = where the oscillator length a ho is a parameter which is fixed by the geometry of the trapping potential and a aa , a ab are the scattering lengths. Such a decay rate should be compared with the creation rate Λ+Λ ab h . However it should be noted that, since the decay rate depends on external parameters as the harmonic trapping frequency and since the scattering length could be controlled by Feshbach resonances, it should be possible to realize a regime where the creation rate is much larger than γ BM and detection of cat states is possible. Let us finally note that by fixing the ratio of 87 Rb -87 Rb interaction to 87 Rb -85 Rb interaction to be 2.13, a parameter accessible in the JILA setup [9], the detection time is twice smaller than the case of a single component BEC.

VI. CONCLUSIONS AND PERSPECTIVES

In this paper we investigated the quantum dynamics of a Bose Josephson junction made of a binary mixture of BECs loaded in a double well potential within the two-mode approximation. We focused on the small tunneling amplitude limit and adopted the angular momentum representation for the Bose-Hubbard dimer Hamiltonian. Perturbation theory up to second order in the tunneling amplitude enabled us perform analytical calculations in the symmetric case where Λ a = Λ b = Λ and K a = K b = K. In this way we obtained the energy eigenvalues and eigenstates, whose knowledge is mandatory in order to investigate the quantum evolution of the number difference of bosons between the two potential wells. In order to study the quantum dynamics more easily and analytically, we restricted to the case in which the two species are equally populated and imposed the condition of equal population imbalance of the species a and b between the two wells. We concentrated on the two following initial conditions: completely localized states and coherent spin states, and found a rich and complex behaviour, ranging from small amplitude oscillations and collapses and revivals to coherent tunneling. Finally, we considered the generation of Schroedinger cat states and pointed out their influence on the momentum distribution through the vanishing of the contrast. We showed that the creation time can be rendered short enough in the case of a mixture in order to render their detection possible. That could be crucial in order to build up an experimental protocol to produce and detect cat states within such systems.

We stress that in this work we have chosen to study the symmetric case. This allowed us to obtain analytical results, while giving rise to the relevant phenomenology which characterizes the physics of the junction. The general case of different couplings between the two bosonic species and/or different populations needs to resort to numerical calculations and will be the subject of a future publication together with a quantitative evaluation of the validity of the two-mode approximation and of the precision of the perturbative results obtained in this work [34]. Another interesting issue which deserves further investigation is a careful analysis of the quantum manifestations of the selftrapping transition and in general of the MQST phenomenon in this more general context.

The complex dynamics of the generalized Bose Josephson junctions investigated in the present paper could be experimentally testable within the current technology. For instance, the JILA group recently [9] succeeded in producing a mixture of 85 Rb and 87 Rb atoms, whose interactions are widely tunable via Feshbach resonances. In particular it is possible to fix the scattering length of 87 Rb as well as the interspecies one and to tune the scattering length of 85 Rb. That allows one to explore the parameter space in a wide range and also to realize the symmetric regime Λ a = Λ b = Λ. Because of the high degree of experimental control, such a setup could be employed to reproduce the phenomenology described in this work.

Appendix: Order k 2 eigenvectors

The eigenvectors of the full Hamiltonian (34), up to order k 2 , are:

h (2) 0,0 = 1 - k 2 λ 2 j 2 (j 2 +
(2) 0,m ± 2 = A m2 0, m ± 2 + B + m2 0, (m 2 + 1) ± + B - m2 0, (m 2 -1) ± + C + m2 0, (m 2 + 2) ± +C - m2 0, (m 2 -2) ± ;

m 2 = 0, 1, 1 2 , ( 66
)
where the coefficients are defined as:

A m2 = 1 - k 2 4λ 2
4j 2 (j 2 + 1) (m 2 ) 2 + j 2 (j 2 + 1) -4 (m 2 ) 4 + 3 (m 2 )

2 4 (m 2 ) 2 -1 2 , (67)
B ± m2 = ± k 2λ j 2 (j 2 + 1) -m 2 (m 2 ± 1) (2m 2 ± 1) , (68)
C ± m2 = k 2 16λ 2

2

  and j b = N b 2 respectively. So the whole Hilbert space has finite dimension, equal to (2j a + 1) ⊗ (2j b + 1) = (N a + 1) ⊗ (N b + 1), thus it depends on the number of bosons of the species a and b respectively. The whole basis {|m a , |m b } is given by the eigenvectors of J a z (J a z |m a = m a |m a ) and J b z (J b z |m b = m b |m b ) with m a = -Na 2 , ..., Na 2 and m

2 FIG. 1 :

 21 FIG. 1: (Color online)Time evolution of the relative boson number difference between the two traps for different timescales. The value of k is k = 0.5 and the boson number is N2 = 10., while λ = 1 on the left panels and λ = 1.3 on the right panels. The case λ = 1. corresponds to the single species BEC, i.e. Λ ab = 0.

2 FIG. 2 : 2 FIG. 3 :

 2223 FIG. 2: (Color online)Time evolution of the relative boson number difference between the two traps for different boson numbers. The value of k is k = 0.1, λ = 1.3 and N2 = 10.(black-dashed line), N2 = 15.(red-dotted line), N2 = 25.(blue-straight line), while θ = π/2 and φ = π/4.

2 FIG. 4 :

 24 FIG. 4: (Color online)Contrast in the time-evolution of J 2x for θ = π/2 and φ = π/4 and even number of bosons. The red line is for N2 = 10., the blue one for N2 = 14. and the black one for N2 = 20.. The interval in which the contrast is zero increases with increasing N2.

j 2 (

 2 Figure 1a (dynamics1_l13.eps)

Acknowledgments

The authors would like to thank M. Salerno for driving their attention on the topic of Bose Josephson junctions and E. Orignac and A. Minguzzi for discussions and for a critical reading of the manuscript.