text:
# Seyfert-Type Dependences of Narrow Emission-Line Ratios and Physical Properties of High-Ionization Nuclear Emission-Line Regions in Seyfert Galaxies

## 1. Introduction

Seyfert nuclei are typical active galactic nuclei (AGNs) in the nearby universe. They have been broadly classified into two types based on the presence or absence of broad (typically $``$ 2000 km s-1) permitted emission lines in their optical spectra (Khachikian, Weedman 1974); Seyfert nuclei with broad lines are type 1 (hereafter S1), while those without broad lines are type 2 (S2). This difference is thought to be due to a dependence of the visibility of broad-line regions (BLRs) on the viewing angle. Accordingly, the BLR is thought to be located in a very inner region (a typical radial distance from the central black hole is $`r0.01`$ pc; see, e.g., Peterson 1993) and surrounded by a geometrically- and optically-thick dusty torus (AGN unified model; see Antonucci 1993 for a review).

Contrary to the BLR emission, narrow (typically $``$ 1000 km s-1) permitted and forbidden emission lines, arising from narrow-line regions (NLRs), are seen in the spectra of both S1s and S2s. Therefore, the NLR is believed to be located far from the nucleus, and thus not to be hidden by dusty tori. If this is the case, the observed physical properties of ionized gas in NLRs do not depend on the viewing angle toward dusty tori. However, several studies have statistically shown that some high-ionization emission lines, such as \[Fe vii\]$`\lambda `$6087, \[Fe x\]$`\lambda `$6374, and \[Ne v\]$`\lambda `$3426, are stronger in spectra of S1s than in those of S2s (e.g., Shuder, Osterbrock 1981; Cohen 1983; Murayama, Taniguchi 1998a; Schmitt 1998; Nagao et al. 2000). Some emission lines whose ionization potential is not so high ($`<`$ 100 eV), but whose critical density is high, such as \[O iii\]$`\lambda `$4363 and \[Ne iii\]$`\lambda `$3869, also show the same Seyfert-type dependence (e.g., Osterbrock et al. 1976; Heckman, Balick 1979; Shuder, Osterbrock 1981; Schmitt 1998; Nagao et al. 2001b). These Seyfert-type dependences in some emission-line intensities seem to conflict with the framework of the AGN unified model. Therefore, some possible models have been proposed to explain such Seyfert-type dependences in the narrow emission-line intensities. One of the proposed ideas is that these Seyfert-type dependences in the strength of the high-ionization emission lines are caused by a viewing-angle dependence of the visibility of highly-ionized dense gas clouds, which are located very close to the nucleus, and thus can be hidden by dusty tori (Murayama, Taniguchi 1998a, b; Nagao et al. 2000, 2001b). On the contrary, Schmitt (1998) pointed out that the Seyfert-type dependences in the intensities of the high-critical-density transition can be understood if the intrinsic difference in the NLR size between S1s and S2s (see, e.g., Schmitt, Kinney 1996) is taken into account.

Which is the origin of the Seyfert-type dependence of those emission lines, the inclination effect (i.e., obscuration by dusty tori) or the intrinsic difference in the size of NLRs? To investigate this issue, it should be examined how emission-line strengths depend on the Seyfert-type. Thus, we have compiled emission-line flux ratios of many Seyfert galaxies from the literature. In this paper, we show the Seyfert-type dependence of various emission-line flux ratios based on this compiled database and discuss the origin of the Seyfert-type dependence.

## 2. The Data

### 2.1. Data

In this section, we briefly describe the properties of the database of emission-line flux ratios. See Nagao (2001) for details of this database. The database contains various emission-line flux ratios of 355 Seyfert nuclei in total; 33 narrow-line S1s (NLS1s; see, e.g., Osterbrock, Pogge 1985), 48 broad-line S1s (BLS1s), 48 Seyfert 1.2 galaxies (S1.2s), 78 Seyfert 1.5 galaxies (S1.5s), 9 Seyfert 1.8 galaxies (S1.8s), 28 Seyfert 1.9 galaxies (S1.9s), 6 S2s with broad emission in near-infrared spectra (S2NIR-BLRs), 12 S2s with broad polarized emission (S2HBLRs), and 93 S2s without any broad emission in their spectra (S2-s). In this paper, S1.2s are basically included in BLS1s, and S1.8s, S1.9s, and S2NIR-BLRs are gathered into the class of “S2 with reddened BLR (S2RBLR)”. When necessary, NLS1s and BLS1s are gathered into the class of “S1total”, and S2RBLR, S2HBLRs, and S2-s are gathered into the class of “S2total”. We adopt the classification of the Seyfert types by Véron-Cetty and Véron (2000).

Since it is often difficult to measure the flux of narrow Balmer components accurately for S1s (and S1.5s), the correction for dust extinction adopting the Balmer decrement method (see, e.g., Osterbrock 1989) would cause possible systematic errors. Therefore, we did not correct any dust extinction effect on the emission-line flux ratios in the database. The effects of dust extinction on the following discussion are mentioned in subsection 3.2.

### 2.2. Selection Biases

The sample used here is not a complete one in any sense. Therefore, it is necessary to check whether or not the sample is appropriate for the statistical and comparative investigations carried out in the following sections. Since some possible biases may be caused if there are large systematic differences in their redshift or intrinsic luminosity distributions, we check their distributions below.

First, we check the frequency distribution of redshift for each type of Seyfert galaxies. The histogram of the frequency distribution, the median, the average, and the 1 $`\sigma `$ deviation of the redshift for each type of Seyfert galaxies are presented in figure 1 and table 1. It appears that the S1s and the S1.5s tend to have higher redshifts than the S2s. In order to check whether or not the frequency distributions of the redshift are statistically different among the various types of Seyfert galaxies, we apply the Kolmogorov–Smirnov (KS) statistical test (see, e.g., Press et al. 1988). The null hypothesis is that the redshift distributions of two types of Seyfert galaxies (in any combination) come from the same underlying population. The derived KS probabilities (i.e., the probabilities of two samples drawn from the same parent population) are summarized in table 2. The KS test leads to the following results: (1) the redshift distributions of the NLS1s, the BLS1s, and the S1.5s are statistically indistinguishable, (2) those of the S2RBLR, S2HBLR, and S2- are also statistically indistinguishable, while (3) the NLS1s, the BLS1s, and the S1.5s have systematically higher redshifts than the S2RBLR, S2HBLR, and S2-. Therefore, we must keep in mind that some properties investigated in the following sections may be affected by the redshift bias. We will check this possibility when we carry out statistical investigations of the emission-line intensity ratios (see subsection 3.2).

Second, we check the intrinsic luminosity distribution of each type of Seyfert galaxies. Following the AGN unified model, the nuclear nonthermal continuum radiation of S2s is absorbed by dusty tori and cannot be observed directly. This results in weaker continuum emission in the sample of S2s than in the sample of S1s, if the distribution of the intrinsic luminosity is the same between the samples of S1s and S2s. In other words, we may pick up intrinsically luminous S2s compared to S1s by survey observations with a certain limiting flux. Since such biases may affect the statistical properties of our samples, we have to check the intrinsic luminosity distribution of each type of Seyfert galaxy using some isotropic radiation. Here we investigate the distributions of the mid- and far-infrared luminosity, i.e., IRAS 25 $`\mu `$m and 60 $`\mu `$m luminosities.111 In this paper, a cosmology of $`H_0`$ = 50 km s-1 Mpc-1 and $`q_0`$ = 0 is assumed. The data of the IRAS 25 $`\mu `$m and 60 $`\mu `$m are taken from Moshir et al. (1992). These luminosities are thought to scale the nuclear continuum luminosity which is absorbed and re-radiated by dusty tori, and to have little viewing-angle dependence (e.g., Pier, Krolik 1992; Efstathiou, Rowan-Robinson 1995; Fadda et al. 1998). Note that the effect of star formation at 25 $`\mu `$m is less than at 60 $`\mu `$m though active starburst may contribute to the 25 $`\mu `$m luminosity, not only the 60 $`\mu `$m luminosity. The histograms of the frequency distribution, the median, the average, and the 1 $`\sigma `$ deviation of the IRAS 25 $`\mu `$m and 60 $`\mu `$m luminosities are given in figures 2 and 3, and table 2. The KS test leads to the result that the samples of S1s and S1.5s have higher IRAS 25 $`\mu `$m luminosity than that of the S2s, though the samples are statistically indistinguishable in the IRAS 60 $`\mu `$m luminosity. Since the intrinsic luminosity of Seyfert nuclei is more accurately represented by the IRAS 25 $`\mu `$m luminosity, this may mean that the samples of S1s and S1.5s are more luminous than that of S2s, statistically. We check whether or not this difference in the intrinsic luminosity distribution affect the later discussions in subsection 3.2.

## 3. Results

### 3.1. Emission-Line Ratios with Balmer lines

Historically, the diagnostic diagrams proposed by Veilleux and Osterbrock (1987; hereafter VO87) have often been used to examine the physical properties of gas clouds in NLRs. The VO87 diagrams are made from the emission-line flux ratios of a forbidden line to a narrow component of a hydrogen Balmer line; e.g., \[O iii\]$`\lambda `$5007/H$`\beta `$, \[N ii\]$`\lambda `$6583/H$`\alpha `$, and so on. However, since it is often difficult to measure the narrow Balmer components for S1s (and S1.5s) accurately, it is unclear whether or not these flux ratios of S1s (and S1.5s) can be used to study the properties of NLRs and can be compared with those of S2s. Therefore, prior to comparing the frequency distributions of these emission-line flux ratios among various types of Seyfert galaxies, we examine whether or not the VO87 diagrams can work even for S1s and S1.5s.

In figure 4, we show the diagnostic diagrams proposed by VO87, i.e., the diagrams of the emission-line flux ratio of \[O iii\]$`\lambda `$5007/H$`\beta `$ versus that of \[N ii\]$`\lambda `$6583/H$`\alpha `$, \[S ii\]$`\lambda \lambda `$6717,6731/H$`\alpha `$, and \[O i\]$`\lambda `$6300/H$`\alpha `$. The compiled data of Seyfert galaxies are plotted in these diagrams with the data of extragalactic H ii systems for references. Note that these emission-line ratios are little influenced by dust extinction, because the wavelength separations between the concerned two lines are small. It appears that the data of the S1s and the S1.5s show larger scatters than those of the S2s in each diagram, especially in the diagram of \[O iii\]$`\lambda `$5007/H$`\beta `$ versus \[S ii\]$`\lambda `$6717,6731/H$`\alpha `$ (figure 4b). Does this suggest that there is a systematic difference in the physical properties of NLR gas clouds between the S2s and the others, or that the measurements of the narrow components of Balmer lines are not well determined for the S1s and the S1.5s? To examine this issue, we investigate the relationship of the emission-line flux ratio between H$`\gamma `$/H$`\beta `$ and H$`\alpha `$/H$`\beta `$ (figure 5). When we assume the case B approximation, their theoretically predicted ratios are H$`\alpha `$/H$`\beta `$ = 2.9 and H$`\gamma `$/H$`\beta `$ = 0.47 (see, e.g., Osterbrock 1989).222 These predictions are insensitive to gas temperature and hydrogen density except for a very high-dense gas clouds as existing in BLRs. In the range of 5000 K $`<`$ $`T`$ $`<`$ 20000 K and 102 cm-3 $`<`$ $`n_\mathrm{H}`$ $`<`$ 104 cm-3, these predicted flux ratios vary within 10 % (see Osterbrock 1989). Note that, in NLRs of Seyfert galaxies, the H$`\alpha `$ emission may be enhanced by collisional excitation (e.g., Ferland, Netzer 1983; Halpern, Steiner 1983). In figure 5, although most of the S2s can be well described by the case B prediction taking the effects of dust extinction into account (approximately 0 mag $`A_V`$ 3 mag), most of the S1s cannot be explained by the case B prediction with dust extinction. This suggests that the fluxes of the narrow components of Balmer lines are not well determined for the S1s (and S1.5s) although other possibilities (e.g., contribution of optically-thin gas clouds) cannot be ruled out. In any case, it is safe to avoid the fluxes of the narrow components of the Balmer lines in order to investigate the physical properties of NLR gas clouds in S1s and in S1.5s and in order to study any systematic differences of the NLRs between S1s and S2s.

### 3.2. Emission-Line Ratios without Balmer Lines

Some other diagnostic diagrams in which the narrow Balmer emission is not used have been proposed to discuss the physical properties of gas clouds in Seyfert nuclei (e.g., Baldwin et al. 1981; Ohyama 1996; Nagao et al. 2001a; see also Halpern, Steiner 1983). Such diagnostic diagrams seem to be highly useful in investigating both the physical properties of NLRs of S1s and S1.5s and any systematic difference from those of S2s. Therefore, we consider statistical properties of the NLRs only by using various forbidden emission-line flux ratios (see also Nagao et al. 2001a).

In figures 6a–p, we show the frequency distributions of the compiled emission-line ratios, \[O i\]$`\lambda `$6300/\[O iii\]$`\lambda `$5007, \[O ii\]$`\lambda `$3727/\[O iii\]$`\lambda `$5007, \[O i\]$`\lambda `$6300/\[O ii\]$`\lambda `$3727, \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007, \[S ii\]$`\lambda `$6717/\[S ii\]$`\lambda `$6731, \[O i\]$`\lambda `$6300/\[S ii\]$`\lambda \lambda `$6717,6731, \[O ii\]$`\lambda `$3727/\[S ii\]$`\lambda \lambda `$6717,6731, \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007, \[O i\]$`\lambda `$6300/\[N ii\]$`\lambda `$6583, \[O ii\]$`\lambda `$3727/\[N ii\]$`\lambda `$6583, \[N ii\]$`\lambda `$6583/\[O iii\]$`\lambda `$5007, \[S ii\]$`\lambda \lambda `$6717,6731/\[N ii\]$`\lambda `$6583, \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O ii\]$`\lambda `$3727, \[Ne v\]$`\lambda `$3426/\[O ii\]$`\lambda `$3727, and \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 for the S1s, the S1.5s and the S2s. In these figures, the distributions of the flux ratios for the NLS1s, the BLS1s, the S2RBLRs, the S2HBLRs, and the S2-s are also shown in the right-hand panels. The frequency distributions of \[O ii\]$`\lambda `$7325/\[O ii\]$`\lambda `$3727, \[S iii\]$`\lambda `$9069/\[S ii\]$`\lambda \lambda `$6717,6731, \[N i\]$`\lambda `$5199/\[N ii\]$`\lambda `$6583, \[N ii\]$`\lambda `$5755/\[N ii\]$`\lambda `$6583, \[Ar iii\]$`\lambda `$7136/\[O iii\]$`\lambda `$5007, \[Fe x\]$`\lambda `$6374/\[O iii\]$`\lambda `$5007, and \[Fe xi\]$`\lambda `$7892/\[O iii\]$`\lambda `$5007 for the S1s, the S1.5s and the S2s are shown in figure 7. Since these line ratios have been measured for a subset of the samples, we do not show their histograms for each subclass of S1s and S2s in this figure. In table 3, the median, the average, and the 1$`\sigma `$ deviation of each emission-line flux ratio for the S1s, the S1.5s and the S2s are given.

In order to check whether or not the frequency distributions of these emission-line flux ratios are statistically different among the S1s, the S1.5s and the S2s, we apply the KS test. The resultant KS probabilities are given in table 4. These results can be summarized as follows. (1) As for the emission-line flux ratios of \[O i\]$`\lambda `$6300/\[O iii\]$`\lambda `$5007, \[O ii\]$`\lambda `$3727/\[O iii\]$`\lambda `$5007, \[O i\]$`\lambda `$6300/\[O ii\]$`\lambda `$3727, \[S ii\]$`\lambda `$6717/\[S ii\]$`\lambda `$6731, \[S iii\]$`\lambda `$9069/\[S ii\]$`\lambda \lambda `$6717,6731, \[O i\]$`\lambda `$6300/\[S ii\]$`\lambda \lambda `$6717,6731, \[O ii\]$`\lambda `$3717/\[S ii\]$`\lambda \lambda `$6717,6731, \[N i\]$`\lambda `$5199/\[N ii\]$`\lambda `$6583, \[O i\]$`\lambda `$6300/\[N ii\]$`\lambda `$6583, \[O ii\]$`\lambda `$3727/\[N ii\]$`\lambda `$6583, \[N ii\]$`\lambda `$6583/\[O iii\]$`\lambda `$5007, \[S ii\]$`\lambda \lambda `$6717,6731/\[N ii\]$`\lambda `$6583 and \[Ar iii\]$`\lambda `$7136/\[O iii\]$`\lambda `$5007, there are little or no systematic differences among the S1s, the S1.5s and the S2s. (2) The emission-line flux ratios of \[O ii\]$`\lambda `$7325/\[O ii\]$`\lambda `$3727, \[N ii\]$`\lambda `$5755/\[N ii\]$`\lambda `$6583 and \[Fe xi\]$`\lambda `$7892/\[O iii\]$`\lambda `$5007 appear to be slightly higher in the S1s than those in the S2s, although these differences are statistically insignificant. (3) It is statistically significant that the emission-line flux ratios of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007, \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O ii\]$`\lambda `$3727, \[Ne v\]$`\lambda `$3426/\[O ii\]$`\lambda `$3727, \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 and \[Fe x\]$`\lambda `$6374/\[O iii\]$`\lambda `$5007 are different between the S2s and the other types of Seyfert galaxies (i.e., the S1s and the S1.5s). Since these 23 emission-line flux ratios are not correlated with the redshift and IRAS 25 $`\mu `$m luminosity as shown in figure 8, the above results seem to be almost free from the redshift and luminosity biases mentioned in subsection 2.2.

As mentioned in subsection 2.1, no reddening correction has been made for all of the compiled emission-line flux ratios. Since it is known that the dust extinction is larger on average in S2s than that in S1s (e.g., Dahari, De Robertis 1988), the above results may be caused by the difference in the amount of extinction between the S1s and the S2s. In order to check whether or not this is the case, we examine the effect of the extinction correction for each emission-line flux ratio adopting the Cardelli’s extinction curve (Cardelli et al. 1989). In table 5, we summarize the correction factors, by which the emission-line flux ratios should be multiplied to be converted into the extinction-corrected values for the case of $`A_V`$ = 1 mag. Note that the mean difference of the amount of dust extinction between S1s and S2s is $``$1 mag (Dahari, De Robertis 1988; see also De Zotti, Gaskell 1985). Since the effect of the extinction correction is too small for the cases of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007 and \[Ne v\]$`\lambda `$3426/\[O ii\]$`\lambda `$3727, the differences in these emission-line flux ratios between the S1s and the S2s cannot be attributed only to the effect of the dust extinction. Furthermore, the differences in the flux ratios of \[Ne iii\]$`\lambda `$3869/\[O ii\]$`\lambda `$3727, \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 and \[Fe x\]$`\lambda `$6374/\[O iii\]$`\lambda `$5007 cannot be also interpreted by the effect of the dust extinction. Therefore, we conclude that the AGN-type dependence of these emission-line flux ratios is due not to the difference in the amounts of extinction, but to some other factors, as discussed later.

Here, we mention that these results are consistent with the previous studies. The excess of the flux ratio of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007 in S1s (and in S1.5s) has been reported by Osterbrock et al. (1976), Heckman and Balick (1979), Shuder and Osterbrock (1981), Cohen (1983), and Nagao et al. (2001b). The excess of the intensities of high-ionization iron emission lines of S1s has been mentioned by Shuder and Osterbrock (1981), Cohen (1983), Murayama and Taniguchi (1998a), and Nagao et al. (2000). Schmitt (1998) noted that S1s exhibit higher ratios of \[Ne iii\]$`\lambda `$3869/\[O ii\]$`\lambda `$3727 and \[Ne v\]$`\lambda `$3426/\[O ii\]$`\lambda `$3727 than S2s.

## 4. Discussion

### 4.1. The Origin of the Seyfert-Type Dependence of the Emission-Line Flux Ratios

The results presented in the last section can be summarized as follows: (1) Most of the emission-line flux ratios which show the Seyfert-type dependence contain a high critical-density ($`10^7`$ cm-3) and/or high ionization-potential ($``$ 100 eV) emission lines, except for the ratio of \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007. (2) On the other hand, the flux ratios which do not exhibit the Seyfert-type dependence do not contain high-ionization emission lines. In table 6, we summarize the ionization potential and the critical density of each emission line used in our analysis. Since the \[Fe xi\]$`\lambda `$7892 emission has a very high ionization potential, its relative intensity could be different between S1s and S2s. However, it is unclear in our sample whether or not the \[Fe xi\]$`\lambda `$7892 emission is stronger in the S1s than in the S2s. One of the reasons for this may be the small number of \[Fe xi\]-detected objects. Further observations will be necessary to confirm the difference in the frequency distribution of the \[Fe xi\]$`\lambda `$7892 intensity between S1s and S2s.

We now consider the origin of the Seyfert-type dependence of those emission-line flux ratios. There are two possible alternatives. One is that the high-ionization emission lines arise mainly from dense gas clouds which are located very close to nuclei (Torus HINER333 The term “HINER” means high-ionization nuclear emission-line region (Binette 1985; Murayama et al. 1998); see Murayama, Taniguchi 1998a, b). Since such a region can be hidden by dusty tori, the visibility of the dense gas clouds may depend on a viewing angle toward the tori. This component may correspond to highly-ionized, dense ($`10^{78}`$ cm-3) gas beside the inner wall of dusty tori (see Pier, Voit 1995). This idea was proposed by Murayama and Taniguchi (1998a) to explain the stronger \[Fe vii\]$`\lambda `$6087 emission in S1s compared to S2s. Taking this Torus-HINER component into account, Murayama and Taniguchi (1998b) showed that the difference in the \[Fe vii\]$`\lambda `$6087 intensity between S1s and S2s can be successfully interpreted by a dual-component photoionization model proposed by them (see also Nagao et al. 2001b). The other idea is that the difference in the emission-line flux ratios reflects the intrinsic difference of NLR properties, e.g., physical size, density, temperature, and so on. Osterbrock (1978) mentioned that the systematic difference in the flux ratio of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007 can be understood assuming $`n_\mathrm{H}10^{67}`$ cm-3 for S1s and $`n_\mathrm{H}<10^5`$ cm-3 for S2s. On the contrary, Heckman and Balick (1979) and Cohen (1983) claimed that the origin of the difference in \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007 between S1s and S2s is attributed not to the density difference, but to the temperature difference; i.e., $`T_\mathrm{e}>2\times 10^4`$ K for S1s while $`T_\mathrm{e}10^4`$ K for S2s. Here, it must be noted that these situations may occur when high-density or high-temperature regions are hidden by any obscuring matter in S2s. Namely, these two scenarios do not necessarily mean the intrinsic difference in the NLR properties. On the other hand, Schmitt (1998) reported that the systematic differences in the ratios of \[Ne iii\]$`\lambda `$3869/\[O ii\]$`\lambda `$3727 and \[Ne v\]$`\lambda `$3426/\[O ii\]$`\lambda `$3727 between S1s and S2s can be explained by taking the physically (i.e., not projected) smaller NLR size of S1s compared to that of S2s, as suggested by Schmitt and Kinney (1996), into account. Since smaller NLRs contain fewer ionization-bounded clouds, which radiate low-ionization emission lines selectively, this results in weaker \[O ii\]$`\lambda `$3727 emission compared to \[Ne iii\]$`\lambda `$3869 and \[Ne v\]$`\lambda `$3426 in S1s (hereafter “smaller NLR model”).

Which scheme is more realistic, the obscuration of the Torus HINER in S2s or the smaller NLR model? In either case, both high-ionization gas clouds and low-ionization ones are necessary to explain the observations, and a certain systematic difference in the relative contribution between the two components can be regarded as the origin of the Seyfert-type dependence of the emission-line flux ratios. However, the two ideas make different predictions. The Torus-HINER model predicts that the intrinsic luminosity of high-ionization lines are brighter in S1s than in S2s because the Torus HINER is assumed to be hidden in the S2s (Murayama, Taniguchi 1998a). On the other hand, the smaller NLR model predicts that the intrinsic luminosities of low-ionization lines are lower in the S1s than in the S2s because the low-ionization lines arise mostly from ionization-bounded gas clouds at outer NLRs. Therefore, the Torus-HINER model predicts similar \[O iii\]$`\lambda `$5007 luminosity among the Seyfert types and higher luminosities of \[O iii\]$`\lambda `$4363, \[Ne iii\]$`\lambda `$3869, and \[Fe vii\]$`\lambda `$6087 in S1s than in S2s, while the “smaller NLR model” predicts lower \[O iii\]$`\lambda `$5007 luminosity in S1s than that in S2s, and similar luminosities of \[O iii\]$`\lambda `$4363, \[Ne iii\]$`\lambda `$3869, and \[Fe vii\]$`\lambda `$6087 among the Seyfert types.

We now discuss observational tests for the two models. Since the frequency distribution of the intrinsic luminosity is different among the S1s, the S1.5s and S2s in our sample, as mentioned in subsection 2.2, these emission-line luminosities should be normalized by IRAS 25 $`\mu `$m luminosity to test the above predictions. In figure 9, we show the frequency distributions of the \[O iii\]$`\lambda `$5007, \[O iii\]$`\lambda `$4363, \[Ne iii\]$`\lambda `$3869, and \[Fe vii\]$`\lambda `$6087 luminosities of the S1s, the S1.5s, and the S2s, which are normalized by the IRAS 25 $`\mu `$m luminosity. The median, the average, and the 1$`\sigma `$ deviation of each relative luminosity for the S1s, the S1.5s, and the S2s are given in table 7. In order to examine whether or not the distribution of the relative strength of the emission lines depends on the Seyfert type, we apply the KS test. The resultant KS probabilities are given in table 8. The KS test leads to the following results:

* The \[O iii\]$`\lambda `$5007 luminosity normalized by IRAS 25 $`\mu `$m luminosity is statistically indistinguishable among the Seyfert types. This is consistent with the Torus-HINER model, but in conflict with the prediction of the “smaller NLR model”.
* The \[O iii\]$`\lambda `$4363 and the \[Fe vii\]$`\lambda `$6087 luminosity normalized by IRAS 25 $`\mu `$m luminosity appear to be higher in the S1s and in the S1.5s than in the S2s, although the statistical significance is low. This agrees with the prediction of the Torus-HINER model, but cannot be understood in terms of the “smaller NLR model”.
* The \[Ne iii\]$`\lambda `$3869 luminosity normalized by IRAS 25 $`\mu `$m luminosity is statistically indistinguishable among the Seyfert types. This seems to disagree with the prediction of the Torus-HINER model.

The first two results seem to support the Torus-HINER model. As for the last result, we have to take into account that the critical density of the \[Ne iii\]$`\lambda `$3869 transition is smaller than those of the \[O iii\]$`\lambda `$4363 and \[Fe vii\]$`\lambda `$6087 transitions. This may be the reason why the frequency distribution of the relative luminosity of the \[Ne iii\]$`\lambda `$3869 emission is indistinguishable among the Seyfert types. We therefore conclude that the Torus-HINER model appears to be much more consistent with the observation than the smaller NLR model. We thus adopt the Torus-HINER model in the following discussion.

One of our interesting results is that the emission-line flux ratio of \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007 exhibits a systematic difference in its frequency distribution between the S1.5s and the S2s. If this difference is also attributed to different contributions of the Torus HINER between the two types, it is suggested that the visibility of the \[O iii\]$`\lambda `$5007 emitting region depends on the viewing angle toward dusty tori. Indeed, the isotropy of the \[O iii\]$`\lambda `$5007 emission has sometimes been called into question, especially in radio-loud AGNs. Jackson and Browne (1990) reported that the \[O iii\]$`\lambda `$5007 luminosity of the narrow-line radio galaxies is lower by 5–10 times than that of the broad-line quasars, matched in redshift and extended radio luminosity. On the other hand, Hes et al. (1993) found that these two kinds of radio-loud AGNs show no difference in the \[O ii\]$`\lambda `$3727 luminosity. These results suggest that the \[O ii\]$`\lambda `$3727 emission is isotropic, but the \[O iii\]$`\lambda `$5007 emission is not (see also Baker, Hunstead 1995; Baker 1997; cf., Simpson 1998). Note that this conclusion is based on an assumption that broad-line quasars and (powerful) narrow-line radio galaxies are intrinsically similar, but different only in the viewing angle (see, e.g., Barthel 1989). Polarized \[O iii\]$`\lambda `$5007 emission has been detected in some radio galaxies (di Serego Alighieri et al. 1997), supporting the scenario that a part of the \[O iii\]$`\lambda `$5007 flux is hidden by the tori. Because such polarized \[O iii\]$`\lambda `$5007 emission has also been detected in a Seyfert 2 galaxy, NGC 4258 (Wilkes et al. 1995; Barth et al. 1999), it is interesting to investigate whether or not the \[O iii\]$`\lambda `$5007 emission has an anisotropic property also in the sample of Seyfert galaxies, although Mulchaey et al. (1994) reported a negative result for this possibility.

If the \[O iii\]$`\lambda `$5007 emission of our sample is affected by an orientation-dependent dust obscuration, the frequency distribution of $`L`$(\[O iii\]$`\lambda `$5007) would show a larger scatter than that of $`L`$(\[O ii\]$`\lambda `$3727) as a result of the obscuration of \[O iii\]$`\lambda `$5007 emission at a large inclination angle (see Kuraszkiewicz et al. 2000). Here, the \[O ii\]$`\lambda `$3727 emission is thought to have no viewing-angle dependence. In figure 9, we show a diagram of $`L`$(\[O iii\]$`\lambda `$5007) versus $`L`$(\[O ii\]$`\lambda `$3727), in which no obscuration effect is found. Recently, Kuraszkiewicz et al. (2000) reported that radio-quiet quasars do not show evidence for an anisotropic property of the \[O iii\]$`\lambda `$5007 emission. Therefore, taking the results of both this study and Kuraszkiewicz et al. (2000) into account, radio-quiet AGNs including quasars and Seyfert galaxies may not have the anisotropic property of the \[O iii\]$`\lambda `$5007 emission, contrary to radio-loud AGNs (see also Mulchaey et al. 1994). Then, why is the frequency distribution of the flux ratio of \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007 different between the S1.5s and the S2s? A possible reason may be contamination of lower-ionization emission-line fluxes arising from circumnuclear star-forming regions associated in S2s, because it is known that circumnuclear star-forming activity tends to associate with S2s more frequently than with S1s (e.g., Thuan 1984; Heckman et al. 1995, 1997). Indeed, there is a weak negative correlation between the flux ratio of \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007 and that of $`\nu F_\nu `$(25 $`\mu `$m)/$`\nu F_\nu `$(60 $`\mu `$m), as shown in figure 11. Since the ratio of $`\nu F_\nu `$(25$`\mu `$m)/$`\nu F_\nu `$(60$`\mu `$m) becomes smaller when the star-formation activity contributes much more to the infrared continuum radiation, this correlation suggests that the star-formation activity may enhance the flux ratio of \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007.

### 4.2. Physical Properties of the Torus HINERs

Although the Torus HINER is regarded as a “highly-ionized and dense” component, there are some remaining issues to be investigated: How dense is the Torus HINER? And, how high is the ionization parameter of ionized gas in the Torus HINER? To examine these issues, we carry out single cloud photoionization model calculations using the spectral synthesis code $`Cloudy`$ version 90.04 (Ferland 1996), which solves the equations of statistical and thermal equilibrium and produces a self-consistent model of the run of temperature as a function of depth into a nebula. Here, we assume that a uniform density, dust-free gas cloud with plane-parallel geometry is ionized by a power-law continuum source. The parameters for the calculations are: (1) the spectral energy distribution (SED) of the input radiation, (2) the chemical composition, (3) the hydrogen density of the cloud ($`n_\mathrm{H}`$), (4) the ionization parameter ($`U`$), and (5) the hydrogen column density. We adopt the following input continuum spectrum: (i) $`\alpha `$ = 2.5 for $`\lambda >10\mu `$m, (ii) $`\alpha `$ = –2.0, –1.5, and –1.0, between 10 $`\mu `$m and 50 keV, and (iii) $`\alpha `$ = –2.0 for $`h\nu >`$ 50 keV, where $`\alpha `$ is a spectral index ($`f_\nu \nu ^\alpha `$). Koski (1978) reported that the optical continua of S2s can be approximated by a stellar contribution diluted by a featureless continuum, with the latter component described by a power law with $`\alpha =1.5\pm 0.5`$ (see also Storchi-Bergmann, Pastoriza 1989, 1990; Kinney et al. 1991). We set the gas-phase elemental abundances to be the solar ones taken from Grevesse and Anders (1989) with extensions by Grevesse and Noels (1993). To see the effects of metallicity, we also calculated for the case of $`Z`$ = 0.5 and 2.0 in addition to $`Z=1.0`$. We performed several model runs covering $`10^2`$ $`\mathrm{cm}^3n_\mathrm{H}10^8`$ $`\mathrm{cm}^3`$ and $`10^{3.5}U10^{1.5}`$ for three kinds of SED. Since an unusually strong \[O i\]$`\lambda `$6300 emission with respect to the \[O III\]$`\lambda `$5007 emission is predicted by the model assuming high-density, ionization-bounded gas clouds, we assumed truncated (i.e., matter-bounded) clouds (see, e.g., Murayama, Taniguchi 1998b). To make the gas clouds matter-bounded ones, the calculations were stopped at a hydrogen column density when a Lyman limit optical depth ($`\tau _{912}`$) reached up to 0.1. In this condition, more than 96 % of hydrogen is ionized at the most outer region of a cloud.

In figure 12, the calculated emission-line flux ratios of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007 and \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 are plotted as a function of $`n_\mathrm{H}`$. In order to explain the observed flux ratio of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007 (= 0.122 $`\pm `$ 0.116 for the S1s), rather high density ($`n_\mathrm{H}10^6`$ cm-3) is necessary, regardless of $`U`$, $`\alpha `$, and the metallicity. To examine the ionization parameter which can consistently explain the observed emission-line flux ratios of \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007 and \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 (0.246 $`\pm `$ 0.145 and 0.122 $`\pm `$ 0.116 for the S1s, respectively) adopting $`n_\mathrm{H}10^6`$ cm-3, we show the calculated flux ratios of \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007 and \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 as a function of $`U`$ for the case of $`n_\mathrm{H}10^6`$ cm-3 in figure 13. This figure suggests that the observed emission-line flux ratios can be described by the models with $`U=10^2`$. This value is almost independent of the input SEDs and metallicities.

However, these estimates ($`n_\mathrm{H}10^6`$ cm-3 and $`U=10^2`$) are based on the assumption that all of the observed emission-line fluxes arise from the Torus HINER. In fact, it is evident that a part of the emission-line flux arises from the outer low-density regions. Since the gas clouds in such regions are expected to emit relatively low-ionization emission-line spectra compared to those in the Torus HINER, the Torus-HINER emission must be diluted by such low-ionization emission-line spectra. Therefore, we should remind that the derived properties of Torus HINERs are lower limits. i.e., $`n_\mathrm{H}>10^6`$ cm-3 and $`U>10^2`$.

These results are not modified significantly if we take a smaller value of $`\tau _{912}`$. In figure 14, we show the calculated emission-line flux ratios of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007, and \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 in the case of $`\tau _{912}=0.01`$ as a function of $`n_\mathrm{H}`$ (only the models with solar abundances are shown). The results of the calculations are nearly the same as those for the case of $`\tau _{912}=0.1`$. On the contrary, relatively smaller flux ratios of \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 are predicted if we take a larger value of $`\tau _{912}`$. For reference we show the calculation results for the case of $`\tau _{912}=1.0`$ in figure 15. Accordingly, larger gas density and/or ionization parameters than those estimated for the case of $`\tau _{912}=0.1`$ are necessary to explain the observations in the cases of $`\tau _{912}>0.1`$. In any case, we conclude that the lower limits of the gas density and the ionization parameter of the Torus-HINER are $`n_\mathrm{H}10^6`$ cm-3 and $`U10^2`$.

## 5. Summary

Based on a compilation of the optical emission-line spectra of Seyfert galaxies, we have investigated the Seyfert-type dependences of various emission-line flux ratios. Our analysis was made using only forbidden emission lines. This method has enabled us to compare the physical properties of the NLR among the various types of Seyfert nuclei (e.g., S1, S1.5, and S2).

In consequence of the statistical comparisons of various forbidden emission-line flux ratios among the Seyfert types, we have obtained the following results:

* As for the emission-line flux ratios of \[O i\]$`\lambda `$6300/\[O iii\]$`\lambda `$5007, \[O ii\]$`\lambda `$3727/\[O iii\]$`\lambda `$5007, \[O i\]$`\lambda `$6300/\[O ii\]$`\lambda `$3727, \[S ii\]$`\lambda `$6717/\[S ii\]$`\lambda `$6731, \[S iii\]$`\lambda `$9069/\[S ii\]$`\lambda \lambda `$6717,6731, \[O i\]$`\lambda `$6300/\[S ii\]$`\lambda \lambda `$6717,6731, \[O ii\]$`\lambda `$3717/\[S ii\]$`\lambda \lambda `$6717,6731, \[N i\]$`\lambda `$5199/\[N ii\]$`\lambda `$6583, \[O i\]$`\lambda `$6300/\[N ii\]$`\lambda `$6583, \[O ii\]$`\lambda `$3727/\[N ii\]$`\lambda `$6583, \[N ii\]$`\lambda `$6583/\[O iii\]$`\lambda `$5007, \[S ii\]$`\lambda \lambda `$6717,6731/\[N ii\]$`\lambda `$6583 and \[Ar iii\]$`\lambda `$7136/\[O iii\]$`\lambda `$5007, there are little or no systematic differences among the S1s, the S1.5s and the S2s.
* On the other hand, it is statistically significant that the emission-line flux ratios of \[O iii\]$`\lambda `$4363/\[O iii\]$`\lambda `$5007, \[S ii\]$`\lambda \lambda `$6717,6731/\[O iii\]$`\lambda `$5007 \[Ne iii\]$`\lambda `$3869/\[O iii\]$`\lambda `$5007, \[Ne iii\]$`\lambda `$3869/\[O ii\]$`\lambda `$3727, \[Ne v\]$`\lambda `$3426/\[O ii\]$`\lambda `$3727, \[Fe vii\]$`\lambda `$6087/\[O iii\]$`\lambda `$5007 and \[Fe x\]$`\lambda `$6374/\[O iii\]$`\lambda `$5007 are systematically higher in the S1s and the S1.5s than in the S2s.

These results can be interpreted as that the flux ratios of a high-ionization and a high-critical-density emission line to a low-ionization emission line are systematically higher in the S1s than in S2s.

There are two ideas which can possibly explain the Seyfert-type dependences of the relative strength of high-ionization emission lines; i.e., the viewing-angle dependence of visibility of Torus HINER and the intrinsic difference of the NLR size between S1s and S2s. To discriminate these two possibilities, we have examined the Seyfert-type dependences of some emission-line luminosities. The results are consistent with the prediction of the Torus-HINER model, but disagree with the “smaller NLR model”. We have concluded that the difference of visibility of Torus HINER between S1s and S2s causes the Seyfert-type dependences of the relative strength of the high-ionization emission lines.

In order to investigate the properties of the Torus HINER, we have compared the relative strengths of high-ionization emission lines with the results of the photoionization model calculations. The estimated lower limit of the density and the ionization parameter are $`n_\mathrm{H}10^6`$ cm-3 and $`U10^2`$. These constraints are almost independent of input SEDs and metallicities of the ionized gas.

We would like to thank Gary Ferland for making his code $`Cloudy`$ available to the public. We also acknowledge the referee, Neil Trentham, for useful comments and suggestions. Alberto Rodríguez-Ardila kindly provided information about the emission-line ratios of some Seyfert galaxies. We thank Martin Gaskell for useful comments. This research has made use of the NED (NASA extragalactic database), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under constructed with the National Aeronautics and Science Administration. This work was financially supported in part by Grant-in-Aids for the Scientific Research (Nos. 10044052 and 10304013) of the Japanese Ministry of Education, Culture, Sports, Science, and Technology.

Note added in proof. — C. M. Gaskell (Ap Letts, 24, 43 ) also pointed out that there is a high-density component whose visibility depends on a viewing angle and which is seen only in S1s. This result was derived by the consideration of the emission-line flux ratio of \[S ii\]$`\lambda `$4074/\[S ii\]$`\lambda \lambda `$6717,6731.


format: text

text:
# N=Z Waiting Point Nucleus 68Se: Collective Rotation and High-K Isomeric States

\[

## Abstract

The structure of collective states and quasi-particle excitations in 68Se is analyzed by using the projected shell model. It is found that the $`g_{9/2}`$ orbitals play a decisive role for the structure of this nucleus. At the prolate minimum, alignment of the low-K states causes the observed first backbending in moment of inertia and we predict a second backbending that has not yet been reported. At the oblate minimum, quasi-particles having high-K can form long-lived isomeric states. The structure of 68Se is discussed with respect to the isomers and their possible impact on the determination of its nuclear mass.

\]

Thermonuclear runaways on the surface of neutron stars have been identified as energy sources for type I X-ray bursts . The explosive conditions resulting from the thermonuclear runaway provide the site for rapid sequences of successive proton captures (rp-process) resulting in the creation of nuclei far beyond 56Ni all the way to the proton rich regions  of the chart of nuclides.

X-ray bursts are observed frequently , yet the nucleosynthesis and the correlated energy generation is not completely understood. In addition to the uncertainty of the astrophysical conditions, network simulations of the rp-process are hindered by the lack of experimental information on the structure of nuclei along the rp-process path . Nuclei of particular interest to the rp-process are the N=Z waiting point nuclei . These are the nuclei where further p-capture typically leads to a p-unbound isotone or where the Q value for the (p,$`\gamma `$) reaction is low making the reverse reaction possible. In order for the rp-process to continue,, it needs to wait for $`\beta `$-decay or for the 2p-capture to bridge the waiting point. The time spent at these N=Z waiting points is strongly dependent on the photodisintegration rates, the $`\beta `$-decay lifetimes, and particularly on the masses of nuclei along the path.

The mass of 68Se is particularly relevant to determining whether further p-capture is likely to occur. There have been several attempts  to measure the half-life of 69Br by fragmentation studies at GANIL and MSU. No 69Br was observed but an upper limit of 150ns was set for the half-life. The implication is that 69Br should be p-unbound. This implies that the rp-process should proceed via 2p-capture or $`\beta `$-decay of 68Se. There are some new mass measurements  of several N=Z nuclei in the A=80 mass region from GANIL. The measurement involves the use of the CSS2 cyclotron at GANIL as a high resolution mass spectrometer. The relative mass differences of unknown and known nuclei are determined from their differences in times of flight in the CSS2. The mass excess measured for 68Se was $`52.347\pm `$0.080 MeV. If one now uses the measured 68Se mass, the implication is that 69Br should be bound by greater than 1 MeV. This result is of course inconsistent with the fragmentation experiments. One possible explanation for the discrepancy in the measurements would be if the GANIL mass measurement was that of an isomeric state in 68Se.

The purpose of this Letter is to investigate the role of the structure of 68Se for the mass measurements and therefore for the rp-process. We analyze the structure of both collective and quasi-particle excitations in 68Se using the projected shell model (PSM)  with the idea of looking for the feasibility of the existence of isomeric states in this nucleus.

Great progress has been made on the theory side by the application of advanced shell model diagonalization methods which can explicitly yield the structure wave functions, and matrix elements for even-even, odd-A, and odd-odd nuclei. One successful example is the application of the spherical pf-shell model developed by the Strasbourg-Madrid group . By employing this model, Langanke and Martínez-Pinedo have demonstrated its strength in understanding the nuclear processes that govern those violent astrophysical phenomena .

The rp-process path involves nuclei which are significantly deformed. At the deformed potential minimum, abnormal high-$`j`$ orbitals such as $`g_{9/2}`$ intrude into the pf-shell near the Fermi level. Therefore, the $`g_{9/2}`$ orbitals have to be included in the model space for the present interest. This strongly suggests proper selection of a shell model basis that is capable of describing the undergoing physics within a manageable space.

The PSM  is a shell model diagonalization method which adopts the bases from the well-established Nilsson model with pairing correlation incorporated into them by a BCS calculation. The selection of the bases is first implemented in the multi-quasiparticle (qp) basis with respect to the deformed BCS vacuum $`|0`$; then the broken rotational symmetry is recovered by angular momentum projection technique  to form a shell model basis in the laboratory frame. Finally a shell model Hamiltonian is diagonalized in the projected space. The problem of dimensionality that is often tied with conventional shell-model calculations does not occur in the PSM.

In the present PSM calculation for 68Se, particles in three major shells ($`N=2,3,4`$) for both neutron and proton are activated. The shell model space includes the 0-, 2- and 4-qp states:

$$\{|0,\alpha _{n_i}^{}\alpha _{n_j}^{}|0,\alpha _{p_i}^{}\alpha _{p_j}^{}|0,\alpha _{n_i}^{}\alpha _{n_j}^{}\alpha _{p_i}^{}\alpha _{p_j}^{}|0\},$$
(1)
where $`\alpha ^{}`$ is the creation operator for a qp and the index $`n`$ ($`p`$) denotes neutron (proton) Nilsson quantum numbers which run over properly selected (low-lying) orbitals. Note that the index $`i`$ and $`j`$ in Eq. (1) are general. For example, a 2-qp state can be of positive parity if both quasiparticles $`i`$ and $`j`$ are from the same major shell; it can also be of negative parity if two quasiparticles are from two neighboring major shells. Positive and negative parity states span the whole configuration space with the corresponding matrix in a block-diagonal form classified by the parity.

As the first step in a PSM calculation, one has to find out where the optimal basis is. The shape coexistence feature in 68Se seems to be firmly suggested by the recent experimental data of Fischer et al. . Two coexisting rotational bands were identified, with the ground state band having properties consistent with collective oblate rotation, and the excited band having characteristics consistent with prolate rotation. Constrained Hartree-Fock calculations by Sarriguren et al.  also found coexistence of two energy minima in 68Se, with the oblate solution as the ground state, in supporting the experimental conclusion. Moreover, the observed bands are found to interact rather weakly indicating a high energy barrier . The theoretical barrier between the two well-separated minima is calculated as about 3 MeV . These evidences suggest a possible simplification in our calculation that we can perform the diagonalization separately at the prolate and oblate minimum by neglecting the coupling between the two minima. Collecting all of the above information, we thus perform the following PSM calculations at $`\epsilon _2=0.28`$ for the prolate states, and $`\epsilon _2=0.24`$ for the oblate states.

It is evident from the experimental data  that the rotational states in 68Se are dominated by quadrupole collectivity and pairing interactions. To describe these properties, it is efficient to use a quadrupole plus pairing Hamiltonian

$$\widehat{H}=\widehat{H}_0\frac{1}{2}\chi \underset{\mu }{}\widehat{Q}_\mu ^{}\widehat{Q}_\mu G_M\widehat{P}^{}\widehat{P}G_Q\underset{\mu }{}\widehat{P}_\mu ^{}\widehat{P}_\mu ,$$
(2)
where $`\widehat{H}_0`$ is the spherical single-particle Hamiltonian which contains a proper spin-orbit force, whose strengths (i.e. the Nilsson parameters) are taken from Ref. . The second term in the Hamiltonian is the quadrupole-quadrupole interaction and the last two terms are the monopole and quadrupole pairing interactions, respectively. The interaction strength $`\chi `$ is determined by the self-consistent relation such that the input quadrupole deformation $`\epsilon _2`$ and the one resulting from the HFB procedure coincide with each other . The monopole pairing strength $`G_M`$ is taken to be $`G_M=\left[18.014.5(NZ)/A\right]/A`$ for neutrons and $`G_M=14.5/A`$ for protons. Finally, the quadrupole pairing strength $`G_Q`$ is assumed to be proportional to $`G_M`$, the proportionality constant being fixed to 0.20 in the present work. These interaction strengths are consistent with the values used in the previous PSM calculations for this mass region .

The eigenvalue equation of the PSM for a given spin $`I`$ takes the form

$$\underset{\kappa ^{}}{}\left\{H_{\kappa \kappa ^{}}^IE^IN_{\kappa \kappa ^{}}^I\right\}F_\kappa ^{}^I=0.$$
(3)
The expectation value of the Hamiltonian with respect to a “rotational band $`\kappa `$$`H_{\kappa \kappa }^I/N_{\kappa \kappa }^I`$ defines a band energy, and when plotted as functions of spin $`I`$, we call it a band diagram . A band diagram displays bands of various configurations before they are mixed by the diagonalization procedure of Eq. (3). Irregularity in a spectrum may appear if a band is crossed by another one(s) at a certain spin.

Fig. 1(a) shows our results for the lowest states of each spin, and at the prolate and oblate minima. They are presented in the plots of moment of inertia vs. square of rotational frequency. Remarkable agreement has been obtained for the oblate band. Our calculations suggest a gradual increase in the moment of inertia extended to higher spin states beyond the current experimental band. The basic features of the moment of inertia in the observed prolate band has also been reproduced. In addition to the sharp backbending at $`I=8`$ that has been observed, we predict a second sharp backbending at $`I=16`$. We note that the highest spin state in the measurement of this band was $`I=14`$. Extension of the current experiment to higher spins can thus test our prediction.

In Fig. 1(b), calculated B(E2) values are plotted for the prolate and oblate bands. These values have not been measured experimentally. We found that in the oblate band, the B(E2) values gradually increase until they are saturated at $`I=14`$. In the prolate band, however, drastic variations are predicted. Sudden drops in the curve appear at spins $`I=8`$ and $`I=16`$, corresponding to the two places where the moment of inertia bends back. From both the moment of inertia and B(E2) values, it is thus interesting to see that the co-existing bands in one nucleus exhibit very distinct characteristics where the prolate band is highly irregular and the oblate one is rather smooth.

These interesting observations can be understood by studying the band diagrams. We found that the $`g_{9/2}`$ orbitals play a decisive role for the near yrast structure of this nucleus. Close to the neutron and proton Fermi levels of 68Se, single qp states of $`g_{9/2}`$ with different K quantum numbers are split by the deformed potential. At the prolate minimum, the Fermi levels are surrounded by the low-K states, whereas at the oblate minimum, they are near the high-K states. The low- and high-K states respond rather differently to the rotation, which is reflected in very different pair alignment processes.

In Fig. 2, various configurations are distinguished by different types of lines. The symbols represent the yrast states obtained after the configuration mixing. There are more than 20 bands in each calculation, but only representative ones are displayed for discussion. Note that for the 2-qp bands, there are two closely-lying bands (a neutron and a proton band) because they nearly coincide with each other for the entire spin region.

Fig. 2(a) shows the band diagram calculated at the oblate minimum. There are two 2-qp bands starting at about 5 MeV. These two 2-qp bands (one neutron and one proton band) are based on the $`g_{9/2}`$ quasiparticles with $`K=7/2`$ and $`9/2`$, coupled to $`K=1`$. As spin increases, they gradually approach the g-band, and finally cross with it between $`I=14`$ and 16. The band crossing is so gentle with a very small crossing angle that one cannot see a band disturbance in the yrast solution. In fact, rather smooth behavior has been predicted for the yrast band in both moment of inertia and B(E2) values (see Fig. 1). On the other hand, the two $`g_{9/2}`$ quasiparticles with $`K=7/2`$ and $`9/2`$ can also couple to a 2-qp state with a total $`K=8`$. The $`K=8`$ bands have similar bandhead energies as the $`K=1`$ bands since they originate in the same two $`g_{9/2}`$ quasiparticles.

Fig. 2(b) presents the band diagram at the prolate minimum. In this case, the two $`g_{9/2}`$ quasiparticles with low-K quantum numbers $`K=1/2`$ and $`3/2`$ lie close to the Fermi energy, and 2-qp bands involving these low-K states behave very differently from the oblate case. It can be seen that the two 2-qp bands (again, one neutron- and one proton-band) based on the low-K $`g_{9/2}`$ quasiparticles (coupled to $`K=1`$) sharply across the g-band at $`I=8`$. Because of this crossing, the g-band is disturbed, resulting in the observed first bandbending (see Fig. 1(a)). The crossing changes also the g-band dominance in the yrast wave functions before the crossing to the 2-qp band dominance after the crossing, which leads to the sudden drop in the B(E2) values (see Fig. 1(b)).

A 4-qp state can be formed by coupling of the neutron and the proton 2-qp states discussed above. As shown by the dashed-dotted curve in Fig. 2(b), the 4-qp band can sharply across the 2-qp bands between $`I=14`$ and 16, and dominants the yrast wave functions thereafter. The crossing will result in the second backbending in moment of inertia, as shown in Fig. 1(a), and a bigger drop in B(E2) values at $`I=16`$, as predicted in Fig. 1(b). The highest spin state of the current experimental data for this band is $`I=14`$ . Extension of a few high spin states will be a crucial test of our predictions, thus the correctness of the band crossing interpretation for 68Se.

The theoretical energy spectrum at the oblate minimum is shown in Fig. 3 in comparison with available data . We predict two nearly-degenerate (only one of them is shown here) high-K bands ($`K=8`$) with a bandhead spin $`I=8`$ and an excitation energy of 5 MeV. They are long-lived isomeric states in the sense that no allowed $`\gamma `$-transition matrix elements of low multipolarity can connect these states to the nearby g-band ($`K=0`$). Our calculation shows the high-K band at a lower excitation energy than any other 2-qp bands. Step-wise transitions via other multi-qp states are also unlikely since there are no other nearby states available. Thus, once the isomeric $`I=8`$ states are populated, the high-K states find no path for further $`\gamma `$-decay.

Our prediction on the high-K isomeric states was based on a realistic calculation that has reasonably reproduced all the known structure data in this nucleus. However, we should mention two possible caveats. First, these isomers are based on 2-qp states of the $`g_{9/2}`$ orbitals. Therefore, correct single-particle energies of the neutron and proton $`g_{9/2}`$ orbitals are important. In the present calculation, we employed the newly-adjusted Nilsson parameters  which were obtained by a best estimation for the proton $`g_{9/2}`$ energy. Thus, further experimental data regarding an accurate proton $`g_{9/2}`$ position are very much desired. Second, neutron-proton pairing has not been explicitly considered in our theory. An interesting aspect in our case is that if the proton-neutron pairing interaction is included in our Eq. (2), the two nearly-degenerate 2-qp states will no longer be proton or neutron 2-qp state, but a superposition of them. The interaction will generally push one of them down, and the other up, thus modifying their positions. Experimentally, if one would observe only one of them, that could be a signature of the proton-neutron pairing effect which may not be sensitively felt by the ground state.

The occurrence of the $`I=8`$ isomers that lie very close to the yrast band may have a significant impact on the mass measurement for this nucleus, and a significant influence on the resulting abundances of the rp-process. The possibility that the GANIL measurement  has hit any isomeric state with an excitation greater than 1 MeV would explain the discrepancy with the fragmentation experiments. It should be noted also that the shape isomer, the lowest state in the prolate minimum, which has not been seen in Ref. , lies possibly at about 1 MeV excitation . Thus, experimentally searching for these isomeric states is very important.

In summary, we calculate the structure of 68Se using the PSM. Our results show excellent agreement with recent experimental findings regarding the spectroscopy of this nucleus: the coexisting oblate and prolate minima, the backbend at $`I=8`$. We further predict an additional backbend at $`I=16`$ and a number of high K isomers at approximately 5 MeV above the ground state. The existence of isomeric states in this nucleus may provide resolution for the observed differences between the different types of measurements, and it may play a significant role in the mass flow and therefore resulting abundances of the rp-process. The latter point will be investigated further via network calculations of the rp-process.

Valuable discussions with J.A. Sheikh and C.J. Lister are acknowledged. Y.S. thanks the colleagues at the University of Notre Dame for their warm hospitality, where this collaboration was initiated. The support of the National Science Foundation under grant PHY 99-01133 and the graduate school of the University of Notre Dame are gratefully acknowledged.


text:
# Mean-field analysis of collapsing and exploding Bose-Einstein condensates

## I Introduction

Bose-Einstein condensates (BECs) of trapped atomic vapor have been realized in several atomic species Anderson ; Davis ; Bradley ; Fried ; Robert ; Modugno . One of the remarkable features of these systems is that both strength Inouye  and sign Cornish ; Roberts01 ; Donley  of interactions between atoms can be tuned by adjusting an external magnetic field near a Feshbach resonance Feshbach . This has opened up the possibility to study collapsing and exploding BECs in a controllable manner, thereby offering a stringent test of the Gross-Pitaevskii mean-field theory GP .

A trapped BEC with attractive interactions may be formed Ruprecht  when quantum pressure arising from Heisenberg’s uncertain principle counterbalances the attractive force between atoms. This is possible when the parameter defined by $`kN_0|a|/d_0`$ is below a critical value $`k_{\mathrm{cr}}`$, where $`N_0`$ is the number of BEC atoms, $`a`$ the s-wave scattering length, and $`d_0`$ the size of the ground-state wave function of a harmonic trap. When the parameter $`k`$ exceeds $`k_{\mathrm{cr}}`$, attractive force dominates quantum pressure, causing BEC to collapse. In the experiments performed by a Rice group Sackett ; Sackett99 ; Gerton01 , BEC atoms are continuously supplied from a supercooled thermal gas, and collapse and growth cycles of BEC have been observed Gerton01 . In the experiments performed at JILA Cornish ; Roberts01 ; Donley , in contrast, $`N_0`$ is fixed and $`a(<0)`$ is decreased by using the Feshbach resonance so as to meet the condition $`k>k_{\mathrm{cr}}`$. This technique enabled them to determine the value of $`k_{\mathrm{cr}}`$ Roberts01 , and to observe exploding atomic ejection from collapsing BEC Cornish ; Donley , a phenomenon called “Bosenova” whose origin is currently under controversy Kagan98 ; SaitoL ; SaitoA ; Duine .

The Gross-Pitaevskii (GP) equation GP  has widely been used to study mean-field properties of BECs, offering fairly good quantitative account of a rich variety of experiments for repulsive BECs Dalfovo . In the case of the attractive BECs, however, it is by no means clear to what extent mean-field theory is valid, for the attractive interaction might enhance many-body quantum correlations. In fact, a measured value of $`k_{\mathrm{cr}}`$ Roberts01  is significantly smaller than that predicted by a mean-field theory Ruprecht ; Gammal , the origin of the discrepancy being not yet understood.

The recent quantitative measurements on collapsing and exploding BEC reported in Ref. Donley  have motivated us to investigate to what extent mean-field theory can explain the experimental observations. This is the main purpose of this paper. By numerically integrating the GP equation with atomic loss, we show that the phenomena reported in Ref. Donley , such as decay of condensates, ejection and refocus of atomic bursts, and jet formation, are reproduced by our numerical simulations. In addition, we predict that various patterns in atomic density are formed in the course of collapse.

This paper is organized as follows. Section II briefly reviews the GP equation with atomic loss due to inelastic collisions and describes the method of analysis of our simulations. Section III reports our results of numerical simulations using the GP equation, and compares them with the experimental data of Ref. Donley . Section IV provides the summary of this paper.

## II The Gross-Pitaevskii equation with atomic loss and methods of analysis

The GP equation describes unitary time evolution of a macroscopic “wave function” $`\psi `$, and conserves the total number of atoms $`|\psi |^2𝑑𝐫`$. In reality, however, atoms are lost from the trap due to the two-body dipolar and three-body recombination losses. These effects may be taken into account by incorporating in the GP equation the imaginary terms describing these inelastic processes Kagan98 :

$`i\mathrm{}{\displaystyle \frac{}{t}}\psi `$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}^2\psi +V_{\mathrm{trap}}(𝐫)\psi +{\displaystyle \frac{4\pi \mathrm{}^2a}{m}}|\psi |^2\psi `$ (1)
$`{\displaystyle \frac{i\mathrm{}}{2}}\left(K_2|\psi |^2+K_3|\psi |^4\right)\psi ,`$
where $`V_{\mathrm{trap}}`$ is the trapping potential, and $`K_2`$ and $`K_3`$ denote two-body dipolar and three-body recombination loss-rate coefficients, respectively. The imaginary terms in Eq. (1) are phenomenologically introduced in order to take account of the fact that the two-body and three-body losses are proportional to the square and cube of the atomic density:

$$\frac{}{t}|\psi |^2𝑑𝐫=(K_2|\psi |^4+K_3|\psi |^6)𝑑𝐫.$$
(2)
We assume that the atoms and molecules produced by inelastic collisions escape from the trap without affecting the condensate. The constants $`K_2`$ and $`K_3`$ include Bose statistical factors $`1/2!`$ and $`1/3!`$, respectively, which are needed for BEC factor , and describe the loss rate per atom.

In the situations we consider the two-body loss can be ignored. The two-body and three-body loss rates are given by $`R_2K_2|\psi |^2`$ and $`R_3K_3|\psi |^4`$, and their ratio by $`R_{23}R_2/R_3=K_2/(K_3|\psi |^2)`$. They have the relation $`R_2R_{23}=K_2^2/K_3`$, which is always $`0.1\mathrm{s}^1`$ for $`{}_{}{}^{85}\mathrm{Rb}`$ Roberts00  and $`{}_{}{}^{7}\mathrm{Li}`$ Gerton . When the density is low and $`R_2R_3`$, $`R_20.1\mathrm{s}^1`$. Then both two-body and three-body losses can be ignored since we shall consider the time scale of $`10`$ ms. When the density is high, the three-body loss dominates the two-body loss, and thus the two-body loss becomes unimportant. In our numerical simulations we shall therefore ignore the two-body loss term in Eq. (1). All results presented below are not affected if the two-body loss is taken into account.

The three-body loss, in contrast, plays a crucial role in determining the collapsing dynamics of BEC. When implosion occurs, the atomic density becomes extremely high, and so does the three-body recombination rate, until it stops the growth of the density. The maximum density in the process of implosion is determined by $`K_3`$ and $`a`$ Kagan98 ; SaitoL ; SaitoA . The values of $`K_3`$ far from the Feshbach resonance have been measured for $`{}_{}{}^{87}\mathrm{Rb}`$ Burt , $`{}_{}{}^{23}\mathrm{Na}`$ Kurn , and $`{}_{}{}^{7}\mathrm{Li}`$ Gerton , and they agree with theoretical predictions Moer ; Fedichev ; Esry  within a factor of ten. Near the Feshbach resonance, however, complicated behaviors of $`K_3`$ are predicted Braaten , with no precise experimental data available Roberts00 .

We performed numerical integration of the GP equation (1) using a finite-difference method with the Crank-Nicholson scheme Ruprecht . Since the peak density changes drastically during implosion SaitoL ; SaitoA , we very carefully controlled the time step to avoid error propagation. Initially we prepared the ground-state wave function by the method in Ref. Edwards  for an initial s-wave scattering length $`a_{\mathrm{init}}`$ and an initial number of BEC atoms $`N_0`$. At $`t=0`$ the interaction is suddenly switched from $`a_{\mathrm{init}}0`$ to $`a_{\mathrm{collapse}}<0`$, inducing collapse of the condensate. We use the parameters of $`{}_{}{}^{85}\mathrm{Rb}`$ and assume the same trap geometry (i.e. radial frequency $`\omega _{}/2\pi =17.5`$ Hz and axial frequency $`\omega _z/2\pi =6.8`$ Hz) as that used in Ref. Donley .

In Ref. Donley , atoms after the collapse are classified into three parts: remnant, burst, and missing atoms. The remnant BEC is a dense atomic cloud peaking around the center of the trap, and the burst is a dilute one that spreads broadly around the remnant BEC. In our simulations, we identify the remnant, burst, and missing atoms as follows. The one-dimensional density distributions defined by $`\rho _{\mathrm{axial}}(z)|\psi |^2𝑑x𝑑y`$ and $`\rho _{\mathrm{radial}}(r)|\psi |^2𝑑z`$ show bimodal structures, one peaking at the center and the other spreading around it (an example is shown in Fig. 5b). We identify the former as the remnant and the latter as the burst, and determine the axial and radial coordinates, $`z_b`$ and $`r_b`$, of the boundaries between the remnant BEC and the burst. We then calculate the number of atoms in the remnant BEC as

$$N_{\mathrm{remnant}}_0^{r_b}2\pi r𝑑r_{z_b}^{z_b}𝑑z|\psi |^2.$$
(3)
The burst atoms are defined by the ones outside the boundary. Because of ambiguities in defining the coordinates of the boundaries, $`N_{\mathrm{remnant}}`$ and $`N_{\mathrm{burst}}`$ defined this way have uncertainties of about $`\pm 0.05N_0`$, while the total number of atoms in the trap $`N_{\mathrm{tot}}N_{\mathrm{remnant}}+N_{\mathrm{burst}}`$ is well defined. We define the number of missing atoms as $`N_{\mathrm{missing}}N_0N_{\mathrm{tot}}`$.

## III Results of numerical simulations

### III.1 Decay of the condensate and local intermittent implosions

We first study the situation in which the scattering length is switched from $`a_{\mathrm{init}}=7a_0`$ to $`a_{\mathrm{collapse}}=30a_0`$ at $`t=0`$, where $`a_0`$ is the Bohr radius. Figure 1 shows time evolutions of the peak height $`|\psi _{\mathrm{peak}}|`$ of the wave function, the total number of atoms in the trap $`N_{\mathrm{tot}}`$, and the number of remnant BEC atoms $`N_{\mathrm{remnant}}`$, where we assume $`K_3=2\times 10^{28}\mathrm{cm}^6/\mathrm{s}`$.

Our numerical simulations show that the condensate first contracts slowly with its peak height $`|\psi _{\mathrm{peak}}|`$ gradually increasing. The total number of atoms remains constant during this process, since the recombination loss is negligible at such low densities. At $`t3.6`$ ms, implosion suddenly occurs in a very localized region (the size of the density spike is $`0.1\mu \mathrm{m}`$ while the size of the BEC cloud is several micrometers) for a very short period of time ($`0.1`$ ms). Furthermore the implosion occurs not just once but many times intermittently for about 10 ms SaitoL ; SaitoA . The three-body recombination loss prominently occurs during the implosion, since the atomic density becomes extremely high. Several tens of atoms are lost in each intermittent implosion, resulting in a step-wise decrease of $`N_{\mathrm{tot}}`$. The decay of $`N_{\mathrm{remnant}}`$ is due to both three-body recombination and atomic burst ejection. We note that the behavior of $`N_{\mathrm{remnant}}`$, shown as the plots in Fig. 1, is very similar to the experimental result (Fig. 1b of Ref. Donley ).

The implosions occur not only at the center of the trap but also at other locations on the trap axis, and more than one density spike is often seen simultaneously. A snapshot of the imploding process with $`a_{\mathrm{init}}=0`$, $`a_{\mathrm{collapse}}=30a_0`$, and $`N_0=15000`$ is displayed in Fig. 2, where the image is taken 2.5 ms after the switch of interactions.

Two large spikes are seen on the trap axis. In an isotropic trap, on the other hand, implosions always occur one by one at the center of the trap.

### III.2 Collapse and decay times

The total number of atoms of the system remains constant for some time after the switch of interactions, and suddenly it begins to decay. We call the time at which the sudden decay begins ‘collapse time’ $`t_{\mathrm{collapse}}`$ (see Fig. 1). Since the collapse time is determined mainly by slow accumulation of atoms towards the center of the trap, $`t_{\mathrm{collapse}}`$ only weakly depends on the value of $`K_3`$. In Fig. 1, $`t_{\mathrm{collapse}}=3.6`$ ms, which agrees with the experimental result of 3.7 ms Donley . Figure 3 shows the dependence of $`t_{\mathrm{collapse}}`$ on $`|a_{\mathrm{collapse}}|/a_0`$, where $`a_{\mathrm{init}}=0`$ and $`N_0=6000`$ Elef .

The plots are in good agreement with the experimental ones (Fig. 2 of Ref. Donley ). For another set of parameters $`a_{\mathrm{init}}=89a_0`$ and $`a_{\mathrm{collapse}}=15a_0`$, we obtain $`t_{\mathrm{collapse}}=15.9`$ ms, which is also consistent with the experimental finding Donley . The collapse time increases for larger values of $`a_{\mathrm{init}}`$, since the atomic cloud spreads more widely and therefore it takes longer time for the cloud to get to the center of the trap.

Fitting

$`N_{\mathrm{remnant}}(t)`$ $`=`$ $`[N_0N_{\mathrm{remnant}}(\mathrm{})]e^{\frac{tt_{\mathrm{collapse}}}{\tau _{\mathrm{decay}}}}`$ (4)
$`+N_{\mathrm{remnant}}(\mathrm{})`$
to the plots in Fig. 1 (dotted curve), we obtain the decay time constant $`\tau _{\mathrm{decay}}3.7`$ ms. For $`a_{\mathrm{init}}=0`$, and $`N_0=6000`$ and $`15000`$, $`\tau _{\mathrm{decay}}`$ is almost the same as the above one, and for $`a_{\mathrm{collapse}}=60a_0`$ it becomes $`\tau _{\mathrm{decay}}3.1`$ ms which agrees reasonably well with the experimental finding of 2.8 ms Donley .

### III.3 Remnant, burst, and missing atoms

Figure 4a shows the fractions of remnant, burst, and missing atoms after the implosion has finished for $`a_{\mathrm{init}}=0`$ and $`a_{\mathrm{collapse}}=30a_0`$.

We note that $`N_{\mathrm{remnant}}`$ is much larger than $`N_{\mathrm{cr}}`$ when $`N_0`$ is large, as observed experimentally Donley . To put it differently, once the condensate expands following the collapse, the critical density will not be reached even when $`N_{\mathrm{remnant}}>N_{\mathrm{cr}}`$. This is because the ratio $`\omega _{}/\omega _z`$ is irrational and therefore the axial and radial refocuses of the burst atoms do not occur simultaneously, with no further implosions occurring.

This result presents a sharp contrast with that for an isotropic case shown in Fig. 4b, where the trap frequency is chosen to be the geometric mean $`(\omega _{}^2\omega _z)^{1/3}`$. The numbers of remnant atoms $`N_{\mathrm{remnant}}`$ are always below $`N_{\mathrm{cr}}`$. This is because all collapsing atoms go to the center of the trap simultaneously, and therefore more implosions occur in an isotropic trap, which increases $`N_{\mathrm{missing}}`$ and decreases $`N_{\mathrm{remnant}}`$. Moreover, the implosions occur also when the burst atoms refocus, since they concentrate at the center of the trap. The data in Fig. 4b are taken before the first refocus. Nevertheless, we note that $`N_{\mathrm{remnant}}`$ is already below $`N_{\mathrm{cr}}`$.

We also note that in the axi-symmetric trap the fractions $`N_{\mathrm{remnant}}/N_0`$, $`N_{\mathrm{burst}}/N_0`$, and $`N_{\mathrm{missing}}/N_0`$ are almost independent of $`N_0`$, particularly for $`N_0>6000`$, which is consistent with the experiments Donley . This is a consequence of the fact that the number of implosions occurring in the collapse is roughly proportional to $`N_0`$ and that the numbers of the burst and missing atoms in each implosion are almost constant.

### III.4 Atomic bursts and “jets”

A burst atom cloud is usually too broadly spread and hence too dilute to be seen. However, at every $`\pi /\omega _{}`$ (or $`\pi /\omega _z`$) the cloud refocuses along the axial (or radial) direction and can be observed. Figure 5a shows the column density seen from the direction perpendicular to the trap axis and Fig. 5b shows the one-dimensional distribution $`\rho _{\mathrm{axial}}(z)|\psi |^2𝑑x𝑑y`$, when the burst atoms focus along the trap axis (corresponding to Fig. 3 of Ref. Donley ).

The s-wave scattering length is switched from $`a_{\mathrm{init}}=0`$ to $`a_{\mathrm{collapse}}=30a_0`$ at $`t=0`$ and the image is taken at $`t=33.6`$ ms $`t_{\mathrm{collapse}}+\pi /\omega _{}`$. There are small peaks in the ridge of focus, which appear to correspond to the shoulders of the central peak in Fig. 3c of Ref. Donley . According to our theory the origin of the atomic burst is a release of kinetic energy of local spikes in the atomic density Kagan98 ; SaitoL ; SaitoA .

As regards the burst energy, it is at present difficult to compare our results with the experimental ones. One reason is that although the burst energy depends not only on $`a_{\mathrm{collapse}}`$ but also on $`K_3`$ SaitoA , experimental values of $`K_3`$ as a function of $`a_{\mathrm{collapse}}`$ are not available. The uncertainty in the boundaries between the remnant BEC and the burst and that in Gaussian fitting of the burst profile (see Fig. 5b) give rise to large errors in numerically determining the burst “temperature”, which also make quantitative comparison between theory and experiment difficult. More work (both theoretical and experimental) needs to be done in order to clarify the situation.

In the experiments Donley , prolonged atomic clouds called “jets” are observed in the radial direction when the collapse is interrupted by switching the interaction from attractive to repulsive. The jets are distinguished from the bursts in that the energy of the former is much lower than the latter and that the direction of the jets is purely radial. In Ref. Donley , the origin of the jets is considered to be the highly anisotropic spikes in the atomic density.

We performed numerical simulations under situations similar to the experimental ones. The s-wave scattering length is switched from $`a_{\mathrm{init}}=0`$ to $`a_{\mathrm{collapse}}=30a_0`$ at $`t=0`$, and changed to $`a_{\mathrm{expand}}=100a_0`$ at $`t=2.5`$ ms (a snapshot before expansion is shown in Fig. 2). Figure 6 shows the gray-scale images of the integrated column density $`\rho (x,z)=|\psi |^2𝑑y`$ after a $`t_{\mathrm{expand}}=3.6`$ ms expansion in (a) and (b), and $`t_{\mathrm{expand}}=5.5`$ ms in (c).

In the image (a), prolonged atomic clouds similar to the jets reported in Fig. 5 of Ref. Donley  appear in the radial direction (the main jet at the center and two tiny jets on either side of it). The images (b) and (c), where the sensitivity of the imaging is increased, show that the jets are interference fringes. The parallel fringe pattern is characteristic of the interference between waves emanating from two point sources Andrews . In fact, there are two spikes in the atomic density that play the role of two point sources of matter waves in Fig. 2. We find that the spacing between the fringes is proportional to $`t_{\mathrm{expand}}`$. This supports our interpretation, since the spacing between the fringes from two point sources is known to be given by $`ht_{\mathrm{expand}}/md`$ Andrews , where $`d`$ is the distance between two point sources. In the experimental images in Ref. Donley  (particularly in Figs. 5d-f), the jets seem to be ejected from the edge of the remnant condensate, where the density is too low to form the spikes, which also supports our interpretation that the jets are the consequence of interference. Thus, the experimental observations of jets indicate that atoms expanding from spikes are coherent, suggesting that the burst atoms are also coherent.

### III.5 Pattern formation

We found in Ref. SaitoL  that various patterns in the atomic density are formed in the collapse processes caused by a sudden switch of interactions from repulsive to attractive, the origin of the pattern formation being attributed to the self-focusing effect of the attractive systems.

Here we predict that a similar pattern formation does occur in $`{}_{}{}^{85}\mathrm{Rb}`$ BEC for experimentally available parameters.

Figure 7 shows pattern formation in axi-symmetric traps, where $`N_0=5\times 10^4`$, and the s-wave scattering length is switched from $`a_{\mathrm{init}}=400a_0`$ to $`a_{\mathrm{collapse}}=310a_0`$. The pancake-shape trap has the ratio $`\omega _z/\omega _{}=\sqrt{8}`$ with the same geometric mean frequencies as the cigar-shape trap. In the cigar-shape trap the cylindrical shell structure is formed, and in the pancake-shape trap the layered structure is formed.

Figure 8 shows time evolution of the integrated column density $`\rho (x,z)=|\psi |^2𝑑y`$ for an isotropic trap, where $`N_0`$, $`a_{\mathrm{init}}`$, and $`a_{\mathrm{collapse}}`$ are the same as in Fig. 7.

We see that the density fluctuations grow to form four concentric spherical shells at $`t=7.7`$ ms, which move inwards and collapse one by one; then at $`t=11`$ ms we see that a new shell is being formed

The resolution of the imaging system in Ref. Donley  (7 $`\mu \mathrm{m}`$ FWHM) is inadequate for observing the patterns in Figs. 7 and 8 (the spacing between the shells is $`1\mu \mathrm{m}`$). Expansion of BEC before imaging will blur out the pattern. In order to observe the pattern formation, therefore, we need to improve the in situ imaging method, or to use larger $`d_0`$ and $`|a|`$ to enlarge the pattern.

## IV Summary

We have studied the dynamics of collapsing and exploding BECs by numerically solving the time-dependent GP equation with atomic loss (1), and compared our results with those of the experiments of Ref. Donley . We find that mean-field theory with atomic loss can account for the following experimental findings: (i) It takes the system a certain time $`t_{\mathrm{collapse}}`$ to undergo a sudden decrease in the number of atoms after the jump of $`a`$. (ii) The number of atoms in the condensate decays exponentially with a decay time constant of a few milliseconds. (iii) The burst atoms are ejected in the collapse process, and refocus after every half trap period. (iv) The fractions of remnant, burst, and missing atoms are almost independent of $`N_0`$, and the number of remnant atoms is much larger than the critical number $`N_{\mathrm{cr}}`$ for large $`N_0`$. (v) The jets are observed when the collapse is interrupted by jumping $`a`$ to a positive value.

We have found that these phenomena are attributed to a rapid sequence of local intermittent implosions, and provided a new interpretation of the jets, i.e., the highly anisotropy of the jets is due to the interference fringes. This suggests that the burst atom cloud is coherent.

The validity of the mean-field GP equation is determined by the gas parameter $`na^3`$, and the depletion is given by $`(na^3)^{1/2}`$. When the implosion occurs, $`n|a|^3`$ becomes $`10^3`$ at the peak density in our simulations, which indicates that the mean-field approximation is still valid at least qualitatively.

Our results presented here suggest that the mean-field approximation can be used to describe the collapsing and exploding dynamics at least qualitatively. A more quantitative comparison between experiments and numerical simulations might reveal effects beyond mean-field approximation. This possibility merits further experimental and theoretical study.

## ACKNOWLEDGMENTS

We thank E. A. Donley for valuable comments. This work was supported by a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, and by the Toray Science Foundation.