format: text

text:
# Degeneracy in Density Functional Theory: Topology in 𝑣- and 𝑛-Space

Present Address: \] Department of Physics, University of Missouri-Rolla, Rolla, MO 65409

## Abstract

This paper clarifies the topology of the mapping between $`v`$\- and $`n`$-space in fermionic systems. Density manifolds corresponding to degeneracies $`g=1`$ and $`g>1`$ are shown to have the same mathematical measure: every density near a $`g`$-ensemble-$`v`$-representable ($`g`$-VR) $`n(𝐫)`$ is also $`g`$-VR (except “boundary densities” of lower measure). The role of symmetry and the connection between $`T=0`$ and $`T=0^+`$ are discussed. A lattice model and the Be-series are used as illustrations.

Density functional theory (DFT) of electrons in an external potential $`v(𝐫)`$ uses the ground-state density $`n(𝐫)`$ as basic variable HK ; KS . $`n(𝐫)`$ uniquely determines $`v(𝐫)`$ (apart from a constant); conversely, for a given number, $`N`$, of electrons, $`v(𝐫)`$ obviously determines $`n(𝐫)`$ uniquely if and only if the ground state is non-degenerate. These facts point to the interesting role of degeneracy in the mapping $`v(𝐫)n(𝐫)`$.

In 1982-3, Levy Levy  and Lieb Lieb  discovered examples of well-behaved density functions, $`n(𝐫)`$, which could not be reproduced as non-degenerate, non-interacting ground-state densities of any $`v(𝐫)`$, and were called non-$`v`$-representable (non-VR). In these examples, the $`n(𝐫)`$ could be reproduced as weighted averages of densities of degenerate ground states corresponding to a $`v(𝐫)`$.

Soon afterwards, Kohn Kohn  adopted a lattice version of the Schrödinger equation, $`𝐫r_l`$ ($`l=1,\mathrm{},M`$), in which $`v(𝐫)`$ and $`n(𝐫)`$ become $`v_l`$ and $`n_l`$ and can be viewed as $`M^{}`$-dimensional ($`M^{}`$-D) vectors: $`\stackrel{}{v}(v_1,\mathrm{},v_M^{})`$, with $`M^{}M1`$, $`v_M=0`$; and $`\stackrel{}{n}(n_1,\mathrm{},n_M^{})`$, with $`n_M=N(n_1+\mathrm{}+n_M^{})`$ footnote . The $`M`$ lattice points are enclosed by boundary points on which $`\stackrel{}{v}=+\mathrm{}`$ and, accordingly, all wave functions vanish. It was shown in Kohn  that in the $`M^{}`$-D $`\stackrel{}{n}`$-space there are finite $`M^{}`$-D regions in which all “points” (densities) are VR.

In 1985, Chayes et al. (CCR) Chayes , using a similar lattice model, proved the following important result: Any well-behaved $`\stackrel{}{n}`$ ($`0<n_l1`$) can be uniquely represented as weighted average of $`g`$ degenerate ground-state densities associated with a $`\stackrel{}{v}`$. We call the special case $`g=1`$ 1-VR, and the general case ($`g1`$) $`g`$-VR. Thus, $`\stackrel{}{n}`$-space is filled by manifolds $`S^g`$, in each of which $`\stackrel{}{n}`$ is $`g`$-VR. Each $`S^g`$ maps on a corresponding manifold $`Q^g`$ in $`\stackrel{}{v}`$-space.

This paper aims at an understanding of the topologies of the regions $`S^g`$ and $`Q^g`$ and of the nature of the boundaries between them. This has mathematical and physical significance, e.g. in the search for self-consistent solutions of Kohn-Sham (KS) equations ullrichkohn . We shall see that $`g`$-VR densities with $`g>1`$ are not mathematically exceptional and, in particular, do not depend on symmetries of the potential. As illustrations we shall present a finite lattice model and the Be-series for $`Z\mathrm{}`$ schipper . We also discuss the relationship between $`T=0`$ ensembles and $`T=0^+`$ thermal ensembles. We limit ourselves to systems that are everywhere nonmagnetic.

Interior regions in $`\stackrel{}{n}`$-space. We first recapitulate Ref. Kohn  for 1-VR densities, and then generalize to $`g`$-VR.

a) $`g=1`$. Since there is a non-degenerate non-interacting ground state with a finite gap, corresponding to some $`\stackrel{}{v}`$, we can use non-degenerate perturbation theory to calculate the first-order density change, $`{}_{}{}^{1}n_{l}^{}`$, due to a weak perturbing potential $`{}_{}{}^{1}v_{l}^{}`$:

$$^1n_l=\underset{l^{}=1}{\overset{M^{}}{}}\chi _{ll^{}}{}_{}{}^{1}v_{l^{}}^{},l=1,\mathrm{},M^{}.$$
(1)
Since, according to HK , the density $`\stackrel{}{n}`$ of a non-degenerate ground state uniquely determines $`\stackrel{}{v}`$ (no arbitrary constant since $`v_M=0`$), the homogeneous set of equations corresponding to (1) has no solution other than $`{}_{}{}^{1}v_{l^{}}^{}0`$, so that (1) can be inverted:

$$^1v_l=\underset{l^{}=1}{\overset{M^{}}{}}\chi _{ll^{}}^1{}_{}{}^{1}n_{l^{}}^{}l=1,\mathrm{},M^{}.$$
(2)
Thus, any first-order change of $`\stackrel{}{n}`$ preserving $`N`$ produces a density which is also 1-VR. This is the case in the entire $`M^{}`$-D neighborhood surrounding the point $`\stackrel{}{n}`$, in which $`\mathrm{\Delta }`$ remains positive.

b) $`g>1`$ wigner . We denote by $`\mathrm{\Psi }_1,\mathrm{},\mathrm{\Psi }_g`$ a set of $`g`$ orthogonal degenerate ground-state wave functions corresponding to $`\stackrel{}{v}^g`$, with common energy $`E`$ and finite gap $`\mathrm{\Delta }`$. The remaining eigenfunctions are $`\mathrm{\Psi }_{g+1},\mathrm{}`$. What conditions must be imposed on infinitesimal \[$`𝒪(\lambda )`$\] potential changes $`{}_{}{}^{1}\stackrel{}{v}`$ which will preserve this degeneracy? In the space of the $`\mathrm{\Psi }_j`$, the perturbed Hamiltonian matrix is

$$H_{ij}=E_i\delta _{ij}+{}_{}{}^{1}V_{ij}^{}i,j=1,2,\mathrm{},$$
(3)
$`E_i=E`$ ($`i=1,\mathrm{},g`$), $`E_iE+\mathrm{\Delta }`$ ($`ig+1`$), and $`{}_{}{}^{1}V_{ij}^{}_{l=1}^M^{}i|{}_{}{}^{1}v_{l}^{}|j`$. Because of the finite $`\mathrm{\Delta }`$, the off-diagonal $`{}_{}{}^{1}V_{ij}^{}`$ with $`ig`$ and $`jg+1`$, or vice versa, can be removed to all orders in $`\lambda `$ by orthogonal transformations leading to the decoupled $`g\times g`$ block Hamiltonian

$$\stackrel{~}{H}_{ij}=E\delta _{ij}+{}_{}{}^{1}V_{ij}^{}+𝒪(\lambda ^2)i,j=1,\mathrm{},g.$$
(4)
Diagonalization by an orthogonal transformation $`T_{ij}=\delta _{ij}+{}_{}{}^{1}t_{ij}^{}+𝒪(\lambda ^2)`$ and the requirement that the eigenvalues remain degenerate gives a transformed Hamiltonian

$$\overline{H}_{ij}\left(T^1\stackrel{~}{H}T\right)_{ij}=(E+{}_{}{}^{1}ϵ)\delta _{ij}+𝒪(\lambda ^2).$$
(5)
On inverting the transformation and noting that every matrix commutes with $`\delta _{ij}`$, one finds $`\stackrel{~}{H}_{ij}=\overline{H}_{ij}+𝒪(\lambda ^2)`$. Comparison with Eq. (4) leads to $`\frac{1}{2}(g1)(g+2)`$ conditions of equal diagonal and vanishing off-diagonal $`{}_{}{}^{1}V_{ij}^{}`$. Thus, every point $`\stackrel{}{v}^g`$ with a $`g`$-fold degenerate ground state and finite $`\mathrm{\Delta }`$ is embedded in a manifold in $`\stackrel{}{v}`$-space of dimension $`D^g=\left[M^{}\frac{1}{2}(g1)(g+2)\right]`$, in which the degeneracy $`g`$ and a finite $`\mathrm{\Delta }`$ are preserved. Ground-state degeneracies in $`\stackrel{}{v}`$-space are “rare” in the above sense.

Each $`\stackrel{}{v}^g`$ gives rise to a set of ensemble densities,

$$\stackrel{}{n}=\underset{j=1}{\overset{g}{}}w_j\stackrel{}{n}_j,0<w_j1,\underset{j=1}{\overset{g}{}}w_j=1,$$
(6)
with $`\stackrel{}{n}_j=\stackrel{}{n}_j\left(R^g\mathrm{\Psi }_j\right)`$, where $`R^g`$ is a $`g`$-D orthogonal transformation \[$`\frac{1}{2}g(g1)`$ parameters\] and the $`w_j`$ are normalized weights ($`g1`$ parameters). Thus, all $`\stackrel{}{n}^{}`$ in a finite $`M^{}`$-D neighborhood enclosing a $`g`$-VR $`\stackrel{}{n}`$, defined by $`\mathrm{\Delta }(\stackrel{}{n}^{})>0`$ and $`w_j>0`$, are also $`g`$-VR: Degeneracy in $`\stackrel{}{n}`$-space is not “rare”.

Boundary surfaces in $`\stackrel{}{n}`$-space. Except for boundaries of $`\mathrm{D}M^{}1`$, $`\stackrel{}{n}`$-space is completely filled by $`M^{}`$-D regions $`S^g`$, with degenerate ground-state levels and a finite positive gap to the nearest excited state.

Apart from “corners” ($`\mathrm{D}M^{}2`$), interior regions $`S^g`$ in $`\stackrel{}{n}`$-space are bounded by $`(M^{}1)`$-D internal and external “surfaces” $`\mathrm{\Sigma }^g`$, see Fig. 1. Each internal $`\mathrm{\Sigma }^g`$ separates two interior regions $`S^g`$ and $`S^{g+1}`$. A point $`\stackrel{}{n}`$ on such a $`\mathrm{\Sigma }^g`$ corresponds to a $`(g+1)`$-fold degenerate ground state but with one of the $`g+1`$ weights equal to zero. On the $`g`$-side of $`\mathrm{\Sigma }^g`$ there is a $`g`$-fold degeneracy and a small gap to the $`(g+1)`$st state opens up from $`0`$. On the $`(g+1)`$-side there is a $`(g+1)`$-fold degeneracy and one of the $`g+1`$ weights starts from $`0`$ and becomes positive.

External boundaries correspond to either one of the $`n_l`$ becoming zero, reflecting particle conservation; or one of the $`n_l`$ becoming $`1`$, reflecting the Pauli principle. For examples of both kinds, see the lattice model below.

Imposed symmetry conditions. Consider a Hamiltonian, invariant under a group $`G`$ of the lattice with irreducible representations $`h`$ with $`\mathrm{D}=g_h`$. Let the ground state have total degeneracy $`g=_hm_hg_h`$ due to $`m_h`$ occurrences of $`h`$. If a small perturbation $`H^{}`$ respecting the group $`G`$ is imposed, how many further conditions must be met by $`H^{}`$ so that the degeneracy is maintained?

We have shown above that the perturbation must have vanishing off-diagonal matrix elements (MEs) between all degenerate ground states, and equal diagonal MEs.

Off-diagonal MEs. The symmetry of $`H^{}`$ assures immediately that all off-diagonal MEs between $`hh^{}`$ vanish tinkham . Off-diagonal MEs between two states of the same $`h`$ involve, by the Wigner-Eckart Theorem, only one independent parameter. The total number of independent parameters characterizing all off-diagonal MEs between the degenerate ground states is therefore $`_hm_h(m_h1)/2`$.

Diagonal MEs. The requirement that all diagonal MEs of $`H^{}`$ are equal is automatically satisfied for states belonging to the same $`h`$, and involves $`_hm_h1`$ conditions for states belonging to different $`h`$.

Thus, to ensure that the $`g`$-fold degenerate ground state maintains its full degeneracy requires the imposition of $`_h[m_h(m_h1)/2+m_h]1`$ conditions.

Next we ask: for a ground state with $`g=_hm_hg_h`$, how many free parameters are there to create ensemble densities invariant under group $`G`$ of the Hamiltonian?

The most general ensemble density obtainable from $`g`$ degenerate ground states $`\mathrm{\Psi }_1,\mathrm{},\mathrm{\Psi }_g`$ is $`_{\mu =1}^gw_\mu n(\mathrm{\Psi }_\mu ^{})`$, $`\mathrm{\Psi }_\mu ^{}=R^g\mathrm{\Psi }_\mu `$. However, the $`\mathrm{\Psi }_\mu ^{}`$ must be chosen in a restricted way to assure that the density has the full symmetry of the Hamiltonian: Let $`\phi _k^{h,\nu _h}`$ ($`k=1,\mathrm{},g_h`$) be a set of orthogonal eigenfunctions, constituting a basis for a given $`h`$ of the group $`G`$, with occurrence $`\nu _h=1,\mathrm{},m_h`$, such that all $`\phi _k^{h,\nu _h}`$ transform in the same way under $`G`$. The most general density of the required symmetry is

$$\stackrel{}{n}=\underset{h,\nu _h}{}w^{h,\nu _h}\underset{k}{}\stackrel{}{n}(R^{\nu _h}\phi _k^{h,\nu _h}).$$
(7)
The total number of free parameters is given by those in the rotation $`R^{\nu _h}`$ plus the number of independent weights $`w^{h,\nu _h}`$, i.e. $`_h[m_h(m_h1)/2+m_h]1`$. We observe that the total number of constraints on $`H^{}`$ is equal to the total number of free parameters in $`\stackrel{}{n}`$. Therefore the dimensionality of the density invariant under $`G`$ associated with any set of degenerate states is the same, independent of $`G`$.

Lattice Example. We illustrate our topological results in an example including symmetry requirements in $`\stackrel{}{v}`$\- and $`\stackrel{}{n}`$-space. Consider the singlet ground state of 4 non-interacting electrons, or, equivalently, 2 spinless non-interacting fermions, whose wave functions $`\phi _l^j`$ are defined on a 13-point fcc lattice (see Fig. 2), and vanish on all exterior points. The discrete Schrödinger equation is

$$\left[\frac{\mathrm{}^2}{2m}_l^2+v_l\right]\phi _l^j=\epsilon ^j\phi _l^j,j,l=1,\mathrm{},13,$$
(8)
where, for fcc lattices, $`_l^2f_l=[_kf_k12f_l]/a^2`$, with $`k`$ running over the nearest neighbors of lattice site $`l`$. Energies and potentials will be measured in units of $`2ma^2/\mathrm{}^2`$.

We consider cases of $`C_{4h}`$ symmetry: $`v_lV_1`$, $`l=\text{1–8}`$; $`v_lV_2`$, $`l=\text{9–12}`$; $`v_{13}=0`$. The lattice wave functions belong to the following representations of $`C_{4h}`$ tinkham : $`A_g`$ ($`s`$-like), $`E_u`$ and $`A_u`$ ($`p`$-like), $`E_g`$ and $`B_g`$ ($`d`$-like), and $`B_u`$ ($`f`$-like). $`A_g`$ occurs three times ($`1A_g`$, $`2A_g`$, $`3A_g`$), $`E_u`$ occurs twice ($`1E_u`$, $`2E_u`$), all other representations occur only once. The eigenfunctions and energies required solutions of 3-, 2- and 1-D secular equations.

Consider the ground state of 2 non-interacting spinless fermions. The lowest level is always $`1A_g`$ ($`1s`$-like). The next level depends on $`V_1`$ and $`V_2`$, see Fig. 3. The $`(V_1,V_2)`$ plane divides into three distinct, infinitely extended regions (shown in yellow), in each of which the eigenfunction of the second level belongs to a particular irreducible representation of $`C_{4h}`$: Region 1a, $`1A_u`$ ($`2p_z`$-like); region 1b, $`2A_g`$ ($`2s`$-like); and region 1c, $`1E_u`$ ($`2p_x`$,$`2p_y`$-like). Inside each region, the ground state is non-degenerate (not counting degeneracies dictated by symmetry, such as in the 2-D representation $`1E_u`$).

On the boundary lines 2a, 2b and 2c (shown in blue), crossings of levels belonging to different representations of $`C_{4h}`$ occur, leading to two-fold (accidental) degeneracies. Along 2a, the potential has cubic symmetry ($`V_1=V_2`$), and the two lowest $`p`$-like levels ($`1A_u`$, $`1E_u`$) coincide. On 2b, the $`2A_g`$ and $`1A_u`$ levels cross, and on 2c, the $`2A_g`$ and $`1E_u`$ levels cross. All three boundary lines meet in a single point, “3”, at $`V_1=V_2=8`$. At this point only, the ground state has a three-fold accidental degeneracy. Note that the degeneracies on 2b and 2c are not due to additional potential symmetries.

As shown above, a ground state with restricted symmetry maintains its degeneracy if $`_h[m_h(m_h1)/2+m_h]1`$ conditions are met. In the present example, $`m_h=1`$ always. Regions 2a, 2b, 2c are therefore lines (1 condition), region 3 is a point (2 conditions) in the $`(V_1,V_2)`$-plane.

We now discuss all $`\stackrel{}{n}`$ with $`C_{4h}`$ symmetry. Like $`\stackrel{}{v}`$, they are fully characterized by two values, $`N_1`$ and $`N_2`$, the total density at points 1–8 and 9–12. $`n_{13}=2N_1N_2`$. $`N_1`$ and $`N_2`$ are confined within a stripe in the $`(N_1,N_2)`$ plane, with an “inner” boundary, with $`n_{13}=1`$, reflecting the Pauli principle ($`N_1+N_2<1`$ would imply $`n_{13}>1`$), and an “outer” boundary, with $`n_{13}=0`$, reflecting particle conservation ($`N_1+N_22`$).

The manifold of allowed densities in the $`(N_1,N_2)`$ plane is subdivided into different regions, corresponding to those of the $`(V_1,V_2)`$-plane. All of the regions in the $`(N_1,N_2)`$ plane have $`\mathrm{D}=2`$. This is in line with our general topological results for $`g`$-VR densities.

Fig. 4 illustrates the highest occupied state along paths in the $`(V_1,V_2)`$ and $`(N_1,N_2)`$ planes (see dotted lines in Fig. 3 top and bottom): $`1A_u`$ with weight 1 (region 1a), $`1E_u`$ with weight 1 (equal-weight combination of $`2p_x,2p_y`$-like states; region 1c), or a linear combination of the degenerate $`1A_u`$ and $`1E_u`$ (region 2a), with fractional weights.

Physical Example: the Be-series. Be has a nominal electronic configuration $`(1s)^2(2s)^2`$ and a term value $`{}_{}{}^{1}S`$. Without interactions, $`2s`$ and $`2p`$ are degenerate. The Be-series consists of neutral Be ($`Z=4`$)rvl , and the 4-electron ions with $`Z=5,6,\mathrm{},\mathrm{}`$. By re-scaling ($`\stackrel{~}{𝐫}=Z𝐫/4`$), this is equivalent to a series of Be-like atoms with $`Z=4`$ but with weakened interaction strength $`\stackrel{~}{e}^2=e^2(4/Z)`$. For $`Z\mathrm{}`$, only the nearly degenerate configurations $`(1s)^2(2s)^2`$ and $`(1s)^2(2p)^2`$, $`(2p)^2=3^{1/2}_{j=1}^3(2p_j)^2`$, need to be considered. This leads to a $`2\times 2`$ secular equation all of whose MEs are exactly known shull . The limiting density is explicitly $`4`$-VR:

$$n_0(\stackrel{~}{𝐫})=2\left[\phi _{1s}(\stackrel{~}{𝐫})^2+0.949\phi _{2s}(\stackrel{~}{𝐫})^2+\frac{0.051}{3}\underset{j=1}{\overset{3}{}}\phi _{2p_j}(\stackrel{~}{𝐫})^2\right],$$
(9)
where the $`\phi `$’s are scaled hydrogenic wave functions (for $`Z=4`$). The configurations $`(1s)^2(2s)^2`$ and $`(1s)^2(2p)^2`$ have amplitudes $`0.974`$ and $`0.225`$. (Unlike in CI expansions bunge , KS densities consist strictly of a finite number of orbital densities for all $`Z`$).

From our previous discussion of boundary surfaces, we observe that the density (9) is not on a boundary, since there is no unoccupied degenerate state. Under these conditions, any density $`n(𝐫)`$ sufficiently near $`n_0(𝐫)`$, in general without any symmetry, is also $`4`$-VR. In particular, there must be a finite range in $`1/Z`$ of the Be-series where densities are $`4`$-VR. This is borne out by numerical calculations savin .

Density ensembles at $`T=0`$ and $`0^+`$. Let $`\stackrel{}{v}^g`$ be a potential with $`g`$ degenerate non-interacting, gapped ground states, $`\mathrm{\Phi }_j`$, generating the $`\frac{1}{2}(g1)(g+2)`$-D set of ensemble densities (6). Now introduce a weak perturbing potential $`{}_{}{}^{1}\stackrel{}{v}`$ of order $`\lambda `$ which, in general, splits the degeneracy and leads to eigenstates and eigenvalues $`\mathrm{\Psi }_j=\overline{R}^g\mathrm{\Phi }_j+𝒪(\lambda )`$; $`E_j=E+{}_{}{}^{1}E_{j}^{}`$, where $`{}_{}{}^{1}E_{j}^{}=𝒪(\lambda )`$, and the corresponding densities are $`\stackrel{}{n}_j`$.

For a given low temperature, $`\beta ^1\mathrm{\Delta }`$, the canonical density corresponding to $`\stackrel{}{v}+{}_{}{}^{1}\stackrel{}{v}`$ is given by

$$\stackrel{}{n}^{}=𝒵^1\underset{j=1}{\overset{g}{}}e^{\beta {}_{}{}^{1}E_{j}^{}}\stackrel{}{n}(\overline{R}^g\mathrm{\Phi }_j)+𝒪(\lambda )+𝒪\left(e^{\beta \mathrm{\Delta }}\right),$$
(10)
where $`𝒵=_j^ge^{\beta {}_{}{}^{1}E_{j}^{}}`$. This density is equal, to leading order, to the density $`\stackrel{}{n}`$ of Eq. (6), if $`{}_{}{}^{1}\stackrel{}{v}`$ is chosen so that $`\overline{R}^g\mathrm{\Phi }_j=\mathrm{\Psi }_j`$ and $`𝒵^1e^{\beta {}_{}{}^{1}E_{j}^{}}=w_j`$. As $`\beta \mathrm{\Delta }\mathrm{}`$, $`{}_{}{}^{1}v\beta ^10`$ and $`\stackrel{}{n}^{}\stackrel{}{n}`$. This establishes the correspondence between ensemble densities at $`T=0`$ and canonical densities at low temperatures, $`\beta \mathrm{\Delta }1`$.

Concluding remarks. The problem of $`v`$-representability Levy ; Lieb  had cast a shadow of uncertainty over DFT, which was very largely clarified in Chayes . Here, building on CCR, we clarify the topology of the $`v(𝐫)`$-$`n(𝐫)`$ mapping. The Be-series and a lattice model are used for illustration. The transition from a dense discrete lattice, $`r_l`$, to a continuous variable, $`𝐫`$, does not appear to offer difficulties. An open issue is extension to spin magnetism.

###### Acknowledgements.

Support by NSF Grants No. DMR-96-30452 and DMR-99-76457 and discussions with A. Savin and C. Umrigar are gratefully acknowledged.


text:
# Constraints from Neutrinoless Double Beta Decay

## Abstract

We examine the constraints from the recent HEIDELBERG-MOSCOW double beta decay experiment. It leads us to the almost degenerate or inverse hierarchy neutrino mass scenario. In this scenario, we obtain possible upper bounds for the Majorana $`CP`$ violating phase in the lepton sector by incorporating the data from the neutrino oscillation, the single beta decay experiments, and from the astrophysical observation. We also predict the neutrino mass that may be measurable in the future beta decay experiments.

The recent neutrino oscillation experiments have shown that neutrinos have masses. On the other hand, the experiments intending to determine directly neutrino mass such as the neutrinoless double beta decay($`(\beta \beta )_{0\nu }`$) and the single beta decay experiments are also on going. In a series of papers, we have discussed the $`CP`$ violation effects in the lepton sector incorporating all these experiments  and the other indirect astrophysical observations . Recently, Klapdor-Kleingrothaus et al. argued the evidence for $`(\beta \beta )_{0\nu }`$ by analyzing the data of the HEIDELBERG-MOSCOW experiment and reported that

$$0.11<m_\nu <0.56\text{eV (95}\%\text{ C.L.)}$$
(1)
with the best fit value, $`m_\nu =0.39`$ eV.

In this paper we reanalyze our studies in response to this announcement, although some papers have discussed the constraints from this data. We especially consider the $`CP`$ violation effects in it. Namely, using the data of Eq.(1), we try to derive the constraints on the $`CP`$ violating phases in the lepton sector by combining the constraints from the neutrino oscillations, the beta decay experiments, and the astrophysical observations. Since the direct test of $`CP`$ violations from a measurement of such as electric dipole moments of leptons seems to be infeasible at present, the $`(\beta \beta )_{0\nu }`$ can be a good channel to detect the $`CP`$ violation effects, although they are indirectly measurable.

The $`(\beta \beta )_{0\nu }`$ experiment gives us the information of the averaged mass $`m_\nu `$ for Majorana neutrinos defined by

$$m_\nu |\underset{j=1}{\overset{3}{}}U_{ej}^2m_j|=|m_1|U_{e1}|^2+m_2|U_{e2}|^2e^{2i\beta }+m_3|U_{e3}|^2e^{2i\rho }|.$$
(2)
Here $`\beta `$ and $`\rho \rho \varphi `$ are $`CP`$ violating phases. The $`U_{aj}`$ is the Maki-Nakagawa-Sakata (MNS) left-handed lepton mixing matrix that combines the weak eigenstate neutrino ($`a=e,\mu `$ and $`\tau `$) with the mass eigenstate neutrino of mass $`m_j`$ ($`j`$=1,2 and 3). The $`U`$ takes the following form in the standard representation :

$$U=\left(\begin{array}{ccc}c_1c_3& s_1c_3e^{i\beta }& s_3e^{i(\rho \varphi )}\\ (s_1c_2c_1s_2s_3e^{i\varphi })e^{i\beta }& c_1c_2s_1s_2s_3e^{i\varphi }& s_2c_3e^{i(\rho \beta )}\\ (s_1s_2c_1c_2s_3e^{i\varphi })e^{i\rho }& (c_1s_2s_1c_2s_3e^{i\varphi })e^{i(\rho \beta )}& c_2c_3\end{array}\right).$$
(3)
Here $`c_j=\mathrm{cos}\theta _j`$, $`s_j=\mathrm{sin}\theta _j`$ ($`\theta _1=\theta _{12},\theta _2=\theta _{23},\theta _3=\theta _{31}`$). Note that three $`CP`$ violating phases, $`\beta `$ , $`\rho `$ and $`\varphi `$ appear in $`U`$ for Majorana particles .

The other experimental constraints on neutrino mass and neutrino mixing angles are as follows: The recent beta decay experiments  restrict another averaged neutrino mass $`\overline{m_\nu }`$ as

$$\overline{m_\nu }\sqrt{\underset{j=1}{\overset{3}{}}|U_{ej}|^2m_j^2}<2.2\text{ eV}.$$
(4)
From the solar neutrino oscillation experiment, we have

$`\mathrm{\Delta }m_{12}^2=m_2^2m_1^2=\mathrm{\Delta }m_{\text{solar}}^2=(220)\times 10^5\text{eV}^2,`$ (5)
$`0.3\mathrm{sin}^22\theta \text{solar}0.93\text{for LMA-MSW},`$ (6)
and

$`\mathrm{\Delta }m_{12}^2=m_2^2m_1^2=\mathrm{\Delta }m_{\text{solar}}^2=(49)\times 10^6\text{eV}^2,`$ (7)
$`8\times 10^4\mathrm{sin}^22\theta _{\text{solar}}8\times 10^3\text{for SMA-MSW},`$ (8)
for the large mixing angle (LMA) and small mixing angle (SMA) MSW solutions, respectively. From the atmospheric neutrino oscillation experiment, we obtain

$$\mathrm{\Delta }m_{23}^2=m_3^2m_2^2=\mathrm{\Delta }m_{\text{atm}}^2=\{\begin{array}{cc}\hfill (17)\times 10^3\text{eV}^2& \text{for normal hierarchy case }\hfill \\ \hfill (17)\times 10^3\text{eV}^2& \text{for inverse hierarchy case}.\hfill \end{array}$$
(9)
The astrophysical observation  gives

$$\underset{i}{}m_i<1.8\text{ eV}$$
(10)
under some reasonable assumptions.

Hereafter we denote the experimental lower and upper bounds in Eq.(1) as $`m_\nu \text{min}`$ and $`m_\nu \text{max}`$, respectively. ($`m_\nu \text{min}m_\nu m_\nu \text{max}`$). Let us first show that the data of the HEIDELBERG-MOSCOW $`(\beta \beta )_{0\nu }`$ experiment in Eq.(1) prefers the almost degenerate or inverse hierarchy neutrino mass scenario for the LMA-MSW solution: Irrespectively of the $`CP`$ violating phases Eq.(2) leads to the inequality that

$$m_\nu m_1|U_{e1}|^2+m_2|U_{e2}|^2+m_3|U_{e3}|^2.$$
(11)
Since we have the constraint $`|U_{e3}|^2<0.03`$ from the oscillation experiments of CHOOZ and SuperKamiokande, Eq.(11) becomes

$$m_\nu <|U_{e1}|^2m_1+|U_{e2}|^2\sqrt{m_1^2+\mathrm{\Delta }m_{12}^2}+0.03\sqrt{m_1^2+\mathrm{\Delta }m_{12}^2+\mathrm{\Delta }m_{23}^2}.$$
(12)
It is apparent from Eqs.(1), (5), (9), and (12) that the normal hierarchy, $`m_1m_2m_3`$, is forbidden. We have no way of distinguishing between the almost degenerate and inverse hierarchy neutrino mass scenarios based on Eq.(1) at this stage. Hence we obtain

$`m_\nu `$ $``$ $`m||U_{e1}|^2+|U_{e2}|^2e^{2i\beta }|,`$ (13)
$`\overline{m_\nu }`$ $``$ $`m,`$ (14)
with $`mm_1m_2`$. Since $`|U_{e3}|^2<0.03`$, $`\mathrm{sin}^22\theta \text{solar}`$ becomes $`4|U_{e2}|^2(1|U_{e2}|^2)`$ and Eq.(13) is rewritten as

$$\mathrm{sin}^2\beta =\frac{1}{\mathrm{sin}^22\theta \text{solar}}\left(1\frac{m_\nu ^2}{m^2}\right).$$
(15)
For LMA-MSW solution, Eq.(15) gives

$$\frac{1}{(\mathrm{sin}^22\theta \text{solar})\text{max}}\left(1\frac{m_\nu ^2}{m^2}\right)\mathrm{sin}^2\beta \frac{1}{(\mathrm{sin}^22\theta \text{solar})\text{min}}\left(1\frac{m_\nu ^2}{m^2}\right).$$
(16)
Here we have denoted the experimental lower and upper limits of Eq.(6) as $`(\mathrm{sin}^22\theta \text{solar})\text{min}`$ and $`(\mathrm{sin}^22\theta \text{solar})\text{max}`$, respectively. The allowed region in the $`\mathrm{sin}^2\beta m`$ plane for the LMA-MSW is shown by the shaded area in Fig.1 with use of the experimental bounds in Eq.(1). Let us superimpose on this allowed region the constraint of the experimental upper bound $`m\text{max}`$ for $`m`$ that is obtained from astrophysical observation and single beta decay. We have $`m<0.6\text{ eV}`$ from Eq.(10) and $`m<2.2\text{ eV}`$ from Eq.(4) using $`\overline{m_\nu }m`$. Namely, at present, the following experimental upper bound is obtained:

$$m<0.6\text{ eV}m\text{max}.$$
(17)
It turns out from Fig.1 that a meaningful bound on the $`CP`$ violating phase $`\beta `$

$$\mathrm{sin}^2\beta \frac{1}{(\mathrm{sin}^22\theta \text{solar})\text{min}}\left(1\frac{m_\nu \text{min}^2}{m\text{max}^2}\right)$$
(18)
is derived for LMA-MSW solution if following condition is satisfied:

$$m\text{max}<\frac{m_\nu \text{min}}{\sqrt{1(\mathrm{sin}^22\theta \text{solar})\text{min}}}.$$
(19)
Thus the constraint on $`\beta `$ from the present experiments is rather weak. So next we show in Fig.2 how this constraint is restricted as the future experiments make progress on the precision measurement, that is, as the lower bounds of $`m_\nu `$ and $`\mathrm{sin}^22\theta _{\text{solar}}`$ are increased. We also obtain a possible lower bound,

$$\frac{1}{(\mathrm{sin}^22\theta \text{solar})\text{max}}\left(1\frac{m_\nu \text{max}^2}{m\text{min}^2}\right)\mathrm{sin}^2\beta ,$$
(20)
if an experimental lower bound $`m\text{min}`$ for $`m`$ (i.e. $`m\text{min}<m`$) is found in the future experiments and the condition $`m_\nu \text{max}<m\text{min}`$ is satisfied. The $`\beta `$ is not restricted in the case of SMA-MSW solution.

Next, following the method used in Ref, we discuss the bound on $`\mathrm{sin}^2\beta `$ by using numerical analysis. In the following discussions, we assume $`m_1m_2m_3`$. The results are scarcely changed for the inverse hierarchical case. In order to obtain the constraints among the observable quantities, let us use $`\overline{m_\nu }`$, $`\mathrm{\Delta }m_{12}^2m_2^2m_1^2`$ and $`\mathrm{\Delta }m_{23}^2m_3^2m_2^2`$ instead of $`m_1`$, $`m_2`$ and $`m_3`$. Namely, Inserting the relations, $`m_2=\sqrt{m_1^2+\mathrm{\Delta }m_{12}^2}`$ and $`m_3=\sqrt{m_2^2+\mathrm{\Delta }m_{23}^2}=\sqrt{m_1^2+\mathrm{\Delta }m_{12}^2+\mathrm{\Delta }m_{23}^2}`$ into Eq.(4) with the unitarity condition that $`|U_{e1}|^2=1|U_{e2}|^2|U_{e3}|^2`$, we obtain the following expressions for $`m_1`$, $`m_2`$, and $`m_3`$ :

$`m_1`$ $`=`$ $`\sqrt{\overline{m_\nu }^2(|U_{e2}|^2+|U_{e3}|^2)\mathrm{\Delta }m_{12}^2|U_{e3}|^2\mathrm{\Delta }m_{23}^2},`$ (21)
$`m_2`$ $`=`$ $`\sqrt{\overline{m_\nu }^2+(1|U_{e2}|^2|U_{e3}|^2)\mathrm{\Delta }m_{12}^2|U_{e3}|^2\mathrm{\Delta }m_{23}^2},`$ (22)
$`m_3`$ $`=`$ $`\sqrt{\overline{m_\nu }^2+(1|U_{e2}|^2|U_{e3}|^2)\mathrm{\Delta }m_{12}^2+(1|U_{e3}|^2)\mathrm{\Delta }m_{23}^2}.`$ (23)
To show $`\beta `$ dependence of $`\overline{m_\nu }`$, we use the center values; $`m_\nu =0.39\text{ eV}`$, $`|U_{e2}|^2=0.29`$ (LMA-MSW solution), $`|U_{e3}|^2=0.03`$, $`\mathrm{\Delta }m_{\text{solar}}^2=4.5\times 10^5eV^2`$, and $`\mathrm{\Delta }m_{\text{atm}}^2=3.2\times 10^3eV^2`$ as a typical case. By inserting Eq.(23) into Eq.(2), we get a relation among $`\overline{m_\nu }`$, $`\beta `$, and $`\rho `$ as depicted in Fig.3. Also from Eq.(23) and Eq.(10), we obtain the upper bound for $`\overline{m_\nu }`$ as

$$\overline{m_\nu }<0.6\text{ eV}.$$
(24)
This has been superimposed on Fig. 3 (b) and (c) giving the upper bound of $`\mathrm{sin}^2\beta `$ as

$$\mathrm{sin}^2\beta <0.7.$$
(25)
Of course, this upper bound depends on the input values of $`m_\nu `$, $`|U_{e2}|^2`$, $`|U_{e3}|^2`$, $`\mathrm{\Delta }m_{\text{solar}}^2`$, and $`\mathrm{\Delta }m_{\text{atm}}^2`$. The Fig.3 also predicts the lower limit of the averaged neutrino mass $`\overline{m_\nu }`$ as $`\overline{m_\nu }>m_\nu `$: Namely, we obtain

$$\overline{m_\nu }>0.39\text{ eV},$$
(26)
indicating that we have a chance for detecting nonzero $`\overline{m_\nu }`$ in the future beta decay experiment.

In conclusion, we have obtained the bounds for the Majorana $`CP`$ violating phases from the recent data of the HEIDELBERG-MOSCOW double beta decay experiment incorporating the data from the neutrino oscillation, the astrophysical observation, and the single beta decay experiments. We have also predicted the lower bound for neutrino mass that may be measurable in the future beta decay experiments.


text:
# Can Reflection from Grains Diagnose the Albedo?

## 1 Introduction

Interstellar dust is poorly understood in spite of its importance within galaxies. The mean extinction law (i.e. scattering plus absorption) for the diffuse Galactic interstellar medium (ISM) is reasonably well known, but the individual absorption and scattering properties are still controversial. Many other basic properties (chemical composition, particle size distribution, particle shape distribution, alignment of the grains in the magnetic field, and others) are also poorly known.

The properties of dust for each scattering event are characterized by two parameters: $`a(\lambda )`$, the albedo or fraction of the extinction that is scattering, and $`g(\lambda )`$, the mean cosine of the scattering angle, averaged over the phase function for single scattering. Geometrical uncertainties prevent the determination of higher moments of the scattering pattern beyond these two. The $`a`$ and $`g`$ are estimated from either reflection nebulae excited by a central star, radiation from a dusty globule reflecting Galactic light, or the “diffuse galactic light” (DGL), the general stellar radiation reflected by the dust distributed throughout the ISM.

Almost always, nebulae have been interpreted with smooth (non-clumpy) models, sometimes with a density gradient. Radiation transfer differs in smooth and clumpy media because photons can escape between the clumps. Uniformity was a reasonable assumption, for three reasons. (a) It is the simplest, and Occam’s Razor has great appeal. (b) Hierarchical clumping had not been convincingly demonstrated by several means. (c) The computational resources required for radiative transport through hierarchical clumps have only recently become widely available.

Any estimate of the albedo is affected by the hierarchical clumping of the ISM. These clumps occur on scales from AU to tens or hundreds of parsecs, as observed in various molecules (CO, formaldehyde, and OH), ionic species, H I, and 100$`\mu `$m (Crovisier & Dickey 1983; Green 1993; Moore & Marscher 1995; Vogelaar & Wakker 1994; Elmegreen & Falgarone 1996; Abergel et al. 1996; Elmegreen 1997; Heiles 1997; Heithausen et al. 1998; Falgarone et al. 1998; Chappell & Scalo 2001; Welty & Fitzpatrick 2001; Faison & Goss 2001; Andrews, Meyer, & Lauroesch 2001). Other galaxies (Stanimirović et al. 1999 for the SMC; Westpfahl et al. 1999 for the M81 group) show the same phenomena in H I. The best observed reflection nebulae are clumpy (e.g. Sellgren, Werner, & Dinerstein 1992; Martini, Sellgren, & DePoy 1999; Knauth et al. 2001). We have examined the images in the SIMBAD database of the reflection nebulae with which we are acquainted (about a dozen); all appear to be clumpy. The observed thermal pressure in the ISM as judged from C I (Jenkins & Tripp 2001) varies by over an order of magnitude within a sample of stars against which the column densities of C I lines can be determined. The hierarchical structure of the ISM is a natural consequence of turbulence, which is scale-free. Turbulent models of the ISM generate hierarchical density structure (e.g., Norman & Ferrera 1996).

Boissé (1990) determined the transport of radiation through a two-phased clumpy medium with isotropic scattering. Witt & Gordon (1996) extended the calculations to anisotropic scattering in a centrally illuminated sphere. They investigated the effects of changing various parameters such as the contrast between the density contrast between the phases or the sizes of the clumps relative to the radius of the sphere. Many of the qualitative effects we find from varying parameters for hierarchical models are discussed by them as well. As the volume fraction of the dense phase was increased, their models varied from widely separated dense clumps up to continuous structures with small holes. Since all cells have either of two densities, we will refer to this type of clumping as “two-phase models”. By contrast, our hierarchical models have an almost continuous distribution of densities.

Witt & Gordon (2000) have investigated radiative transfer in galaxies by means of a specific two-phased ISM, one being 100 times as dense as the other and occupying 15% of the volume. In order to compare two-phase models with hierarchical, we adopt this recipe, with 16 cells along the radius of the sphere (Witt & Gordon used 15). All of our two-phase models assumed ($`a,g)`$ = (0.6, 0.6), which are values typical for optical wavelengths suggested by observations (see plots in Witt & Gordon 2000).

The fact that the ISM is hierarchically clumped is important for the propagation of scattered light, since such a structure has relatively open spaces through which the radiation can move rather freely. The main thrust of this paper is not to try to put limits on the grain properties of real reflection nebulae, but to show the variation among hierarchical models as they are viewed from various angles. We will not be very concerned with the averaged properties of the models, since each real object is viewed from only one direction.

## 2 Hierarchically Clumped Models

The density structure of the ISM is often described as “fractal”, meaning self-similar on all size scales. We will consider models that are hierarchically clumped instead, meaning that they are self-similar over a limited range (about a factor of ten) in sizes. We use a procedure similar to that in Elmegreen (1997): (a) Consider a “supercube”, a portion of which will represent a spherical reflection nebula. The supercube of size L on a side consists of 64 cubical cells stacked along each dimension. For each cell we determine the local density of dust, as explained below. (b) Place $`N`$ points randomly within the supercube. We used $`N`$ = 32. (c) Randomly cast another $`N`$ points, all within a distance $`L/(2\mathrm{\Delta })`$ in each Cartesian axis from each of the points cast in the preceding round. Here the distance $`\mathrm{\Delta }`$ is related to the “fractal dimension”, $`D`$, by $`D\mathrm{log}(N)/\mathrm{log}(\mathrm{\Delta })`$. Allow any of the points that fall outside of the supercube to remain there. (d) Repeat procedure (c), above, twice more, so that the total number of points cast is $`N^4`$. (e) Shift the points outside of the supercube to within it by translating each Cartesian coordinate outside of the supercube by $`L`$ until it lies within. This procedure reflects the points so that if they originally fall off the left side of the supercube they reappear the same distance from the right within the supercube, and correspondingly for up/down and front/back. The density within the supercube is then proportional to the number of points within each cell. (f) Inscribe a sphere within the supercube and place a point source of radiation at its center. The constant of proportionality between the number of points in each cell and the optical depth within the cell is chosen to make the radial optical depth, averaged over all directions, be $`\tau _0`$.

The main parameter of this procedure is $`D`$, which describes the degree of clumping. We considered $`D`$ = 2.3 and 2.6, in the range observed in actual clouds. Another observable parameter is $`\beta `$, the exponent of the power of the projection of the density within the supercube onto the plane of the sky444 The structure function $`F(\delta 𝐫)`$ is defined by $`F(\delta 𝐫)(\mathrm{area})^1I(𝐫)I(𝐫+\delta 𝐫)𝑑x𝑑y`$, and $`P(k)`$ is its Fourier transform, with $`k=|\delta 𝐫|^1`$. Here $`𝐫`$ and $`\delta `$r are vectors, and ($`x,y`$) Cartesian coordinates, in the plane of the sky., so that $`P(k)k^\beta `$, where $`k`$ is the wavenumber. Observations (e.g. Stanimirović et al. 1999) show $`\beta 3`$. With our recipe we find $`\beta =2.8`$ with $`D`$ = 2.6 and $`\beta 2.3`$ with $`D`$ = 2.3. However, none of the conclusions in this paper depend upon which set of hierarchical models we consider. Many different models can have the same $`\beta `$. No one parameter, either $`D`$ or $`\beta `$, can completely describe even the projection of the density distribution of the ISM onto the sky.

The above procedure with the parameters we adopted results in clouds that are quite full of holes, with $``$15% of the projected density distribution being almost blank. The large-scale ISM as shown by H I is similar if a low column density of H I ($`<5\times 10^{19}`$ H atoms cm-2) is taken to be “blank” (i.e., a column density much smaller than the average). Elmegreen (1997) has discussed the emptiness of the ISM, with the large-scale structure in mind. Images of actual reflection nebulae (see Sellgren, Werner, & Dinerstein 1992 for NGC 2023 and NGC 7023) show strong filaments and structure, but there is material at all points within the nebula. For this reason, we sometimes add a uniform density to all cells in the supercube. This constant density does not change the value of $`\beta `$, since the Fourier transform of a constant is a Dirac $`\delta `$-function at the origin, $`k=0`$. Choosing a constant density provides a minimum projected density.

The detailed placement of the clumps of dust also influences the results of our radiative transfer calculations. The clumping is mainly formed by the first round of random casting of 32 points that are subsequently spread by three more rounds of casting points in their vicinities. The locations of the first round of points is uniquely determined by the integer seed of the random number generator555We used ran2(iseed) as described in Press et al. (1992)., and the locations of all subsequent points follow from the initial value of this integer in a complicated but unique way. The response of the model nebula to the central star is strongly affected by the precise placement of the dust relative to the star and, thereby, to the value of the initial seed. We considered suites of models that differed only in their initial seeds. Our procedure amounts to assuming that the star does not affect the density distribution of the dust in its immediate vicinity, so the initial seed is allowed to control the placement of the dust both near and far from the star. To test the sensitivity of results to this assumption, we also considered models with a cavity in the dust distribution within the inner 10% of the radius of the sphere.

The radiative transfer was performed by the Monte Carlo code described by Wood & Reynolds (1999). It involves considering photon packages propagated in a random direction from the central star. We assumed $`a`$ and $`g`$, along with the scattering phase function of Henyey & Greenstein (1941). In most cases we used $`5\times 10^6`$ photons for each model, after checking for a few cases that the results were the same as from a model using $`2\times 10^7`$ photons. Unlike the situation for placing the dust, the initial seed for the radiative transfer makes no difference, as expected from such a large number of stellar photons propagated in random directions.

We calculated the fluxes of scattered radiation and also of starlight, both of which depend on the angle from which the sphere is viewed because of the strong randomly placed clumping. Both fluxes are expressed relative to the stellar flux that would be observed from the unobscured star. We express the fluxes as extinction and scattering optical depths, $`\tau _{\mathrm{ext}}`$ $`\mathrm{ln}(F_{}/F_0)`$ and $`\tau _{\mathrm{sca}}`$ $`\mathrm{ln}(F_{\mathrm{sca}}/F_0)`$, where $`F_0`$ is the flux the star would have if there were no reflection nebula. Each flux or optical depth depends on the direction from which the sphere is viewed. We considered 18 evenly spaced values of cos($`\theta `$) and 36 values of $`\varphi `$, where $`\theta `$ is the polar angle of the sphere and $`\varphi `$ the azimuthal. $`F_{\mathrm{sca}}`$ and $`F_{}`$ can be determined directly (with difficulty!) in real nebulae. $`F_0`$ can be determined if there is enough wavelength coverage to determine the extinction from extrapolation to large wavelengths, from colors if the reddening law is assumed, or if the total FIR flux from the nebula is measured. Otherwise, the ratio of the stellar to the nebular flux gives the difference of $`\tau _{\mathrm{ext}}`$ and $`\tau _{\mathrm{sca}}`$. We now ask, how well can we determine ($`a,g`$) from an observation of ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$)?

## 3 Results

As regards radiative transfer, hierarchical clumping introduces major differences from a uniform distribution or, to a lesser extent, from two-phase models. The main properties of hierarchical models are:

(a) When viewed from various directions, hierarchically clumped nebulae exhibit $`\tau _{\mathrm{ext}}`$ and $`\tau _{\mathrm{sca}}`$ spanning a wide range. A real reflection nebula represents viewing a collection of dust from one particular direction. We can calculate mean properties of models as averages over all viewing directions, but mean properties are relevant only if a large collection of nebulae (such as in another galaxy) are observed. Angle-averaged properties of any model have very limited applicability for interpreting a given reflection nebula. The two-phase models models of Witt & Gordon (1996) also show a wide range of $`\tau _{\mathrm{ext}}`$ at a given $`\tau _0`$.

(b) As expected, the global parameters of the models, such as $`\tau _0`$ and $`D`$, have a significant influence on the radiation transfer. We considered both purely hierarchical models and some with 33% of the mass in a constant density sphere upon which the clumps are superimposed. This uniform background density can have an important effect on the radiation transfer in optically thick nebulae, when the diffusion depends more upon the transparent regions than on the opaque. This rule is familiar in connection with the use of the Rosseland mean opacity in the optically thick limit of radiation diffusion in stars. Whether such a uniform layer exists in reflection nebulae is not clear. The regions between clumps on large (molecular cloud) scales have a vastly lower density than the mean (e.g., Elmegreen 1997).

(c) Models with the same global parameters ($`\tau _0`$, $`D`$, $`a`$, $`g`$), but differing in the placement of the dense clumps relative to the central star, show strong variations in even the averaged fluxes. There are significant differences between hierarchical models in which the star happens to lie within a dense clump of dust, in contrast to those in which it falls within a void (or if there is a central hole). If the star happens to lie within a void, there are lightly reddened paths that reach from the star to the edge of the sphere. If the star is almost unattenuated as seen from a particular direction, the nebula would be interpreted as having the star in front of the dust, except that there can be a significant amount of scattered light that would not be there if the star were truly foreground. In these cases, the albedo derived from uniform models would be $`>1`$.

The power law index, $`\beta `$, of the projected density distribution of our two-phase models is $``$ 0, reflecting the fact that their projected density is non-hierarchical. The two-phase models always provide a smaller range in $`\tau _{\mathrm{sca}}`$($`\theta ,\varphi `$) (about $`\pm `$0.12) than the hierarchical, but their spread in $`\tau _{\mathrm{ext}}`$ is comparable to hierarchical models at low $`\tau _0`$ ($`2`$). The spread at large $`\tau _0`$ decreases because the interclump medium becomes opaque,so that in the ultraviolet (UV) two-phase models make a prediction of $`a`$ that is rather independent of viewing angle.

Figure 1 shows $`\tau _{\mathrm{sca}}`$ plotted against $`\tau _{\mathrm{ext}}`$, for the case $`\tau _0`$ = 2, $`D`$ = 2.6, $`a`$ = 0.6, $`g`$ = 0.6 for a hierarchical distribution of clumps with an initial seed for the density distribution that provides a typical (defined below) placement with respect to the star. Each point represents a particular viewing direction. One sees a large spread in the values of $`\tau _{\mathrm{ext}}`$. The values for $`\tau _{\mathrm{ext}}`$ for this model range up to 6.8, and $`\tau _{\mathrm{sca}}`$ to 2.4. Thus, this one model of clumpy dust can produce a wide range of stellar fluxes and scattered light.

We represent the direction-averaged escaping flux (not mean optical depth) by $`\tau _{\mathrm{ext}}`$, defined by

$$\mathrm{exp}(\tau _{\mathrm{ext}})=\underset{\theta ,\varphi }{}\mathrm{exp}[\tau _{\mathrm{ext}}(\theta ,\varphi )]/N_{\theta ,\varphi },$$
(1)
where $`N_{\theta ,\varphi }`$ is the number of bins we are using ($`=18\times 36`$), all of equal solid angle. If the star is embedded within a clump, the uniform dust component is of reduced importance because scattering occurs within the clump. Similarly, $`\tau _{\mathrm{sca}}`$ represents the angle-averaged scattered radiation, with $`\tau _{\mathrm{sca}}`$ replacing $`\tau _{\mathrm{ext}}`$ in the above relation. The filled square in Figure 1 shows $`\tau _{\mathrm{ext}}`$ and $`\tau _{\mathrm{sca}}`$, which are useful for contrasting the results of models with the same value of $`\tau _0`$ but different initial seeds. As the figure shows, the ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) point is of very limited value in guessing what actual optical depths a given nebula might show if it were viewed from various angles. In this model $`\tau _{\mathrm{ext}}`$, which describes the average optical depth of the stellar flux that escapes from the nebula, is only 0.5$`\tau _0`$, the averaged optical depth through both the clumps and interclump material. It is sobering to recall that only a handful of reflection nebulae have been well observed, while each point in the figure represents a view of a single model with a particular clumping and $`a`$, $`g`$, and $`\tau _0`$.

The solid lines in Figure 1 are the tracks of a spatially uniform dust distribution for the same $`g`$ (= 0.6) and various values of $`a`$, as labeled. These lines were computed with the surprisingly accurate analytical approximation (the scattered flux is underestimated by $``$3% for $`\tau _{\mathrm{ext}}`$ $`<`$ 4) given by Code (1973). Most reflection nebulae have been interpreted using these uniform models, in which case $`\tau _{\mathrm{ext}}`$ is the true optical depth of the nebula. For our hierarchical model, in about 15% of viewing directions we would conclude that the star is in front of most of the dust because the uniform models would give $`a>1`$. In 40% of the directions we would conclude that $`a>0.6`$, the value assumed in the hierarchical model.

Two UV observations of NGC 7023 are shown in the figure. The central star, HD 200775, has foreground extinction from ordinary ISM in addition to the extinction provided by the nebula (see Witt et al. 1993 for discussion and other references to observations of NGC 7023). The circled asterisk marks the nebular parameters derived from Voyager-2 and Hopkins Ultraviolet Telescope (HUT) data if 40% of the total extinction arises within the nebula, and the circled dot if the nebula provides 50%. The dashed lines will be discussed in §4 in connection with the UV DGL.

Figure 2 is the same model as Figure 1, except 33% of the dust is in uniform distribution upon which the hierarchical clumping is superimposed. The axes of the figures are the same as Figure 1 to emphasize that a uniform density increases the scattered light, so long as its optical depth along a radius is not too large ($`1)`$, and imposes a minimum value of $`\tau _{\mathrm{ext}}`$. The range of $`\tau _{\mathrm{sca}}`$ for a given $`\tau _{\mathrm{ext}}`$ is smaller than in Figure 1, but this is an artifact of the uniformity of the extra component that is the only difference in the two figures.

Wide areas of the ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) plane are covered by each hierarchical model, so a given ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) observation can be covered by a wide variety of models. By contrast, uniform models have very limited dependence on $`g`$, and so contain only one effective parameter, $`a`$. The ambiguity of hierarchical models is suggested by Figure 3. It has the same axes as Figures 1 and 2, but only outlines around the areas covered by the points of five hierarchical models. The solid line has a high $`a`$ (= 0.8), moderate $`g`$ (0.6), and low $`\tau _0`$ (0.5). The maximum values of $`\tau _{\mathrm{ext}}`$ are comparatively small because of the low $`\tau _0`$. The dashed line encloses the points for the same model with 33% of the dust mass in a constant density; the scattered points are lower (there is more scattering) because there are no almost empty paths through the nebula. The dot-dashed and dotted lines are the boundaries of the points in Figures 1 and 2, which are for ($`a,g,D,\tau _0)`$ = (0.6, 0.6, 2.6, 2). The long dashed lines enclose the very wide boundaries for ($`a,g,D,\tau _0)`$ = (0.8, 0.85, 2.3, 4). We can have large $`g`$ or moderate, large $`a`$ or moderate, and large $`\tau _0`$ or small and still fit the observations. However, we will see that a minimum $`a`$ of $``$0.5 is required to produce enough scattered light.

Figure 4 shows the averaged ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) values for all of the 21 different initial seeds that we tried. Figure 4$`a`$ shows purely hierarchical models; 4$`b`$, models with 33% of the dust in a constant distribution. The parameters are the same as in Figures 1 and 2: $`a`$ = $`g`$ = 0.6, $`\tau _0`$ = 2. The NGC 7023 points are as before. The open squares are models with $`D`$ = 2.3 instead of 2.6. Within each panel, the $`D`$ = 2.3 models are higher because the radiation, both scattered and stellar, can escape more easily from the more strongly clumped structure. The lines show uniform models, with albedos marked in 4$`a`$. The uniform model appropriate to the actual albedo assumed in the models ($`a`$ = 0.6) is just below the bottom of the figure. The dashed line is the uniform model with $`a`$ = 0.5 and $`g`$ = 0.4 instead of 0.6. The differences are not large in comparison to the effects of the other parameters. The effects on hierarchical models of changing $`g`$ are similar.

Perhaps the most striking difference between the two panels is the lower values of $`\tau _{\mathrm{sca}}`$ for hierarchical models with 33% of the dust in a uniform component. The increase of scattering from the dust between the clumps causes the decrease in $`\tau _{\mathrm{sca}}`$ and greatly reduces the differences between models with different spatial distributions of dust. In either panel, the points from various viewing angles of the two values of $`D`$ are completely intertwined. Reflection nebulae are poor diagnostics of $`D`$ as well as other properties of the ISM.

The three open circles in Figure 4 are two-phase models with three values of the initial seed. Two of these values show the extrema in ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) among the 21 initial seeds that we tested. The tightness of the mean optical depths of the two-phase models models shows the importance of hierarchical geometry, as opposed to simple two-phase clumps. The $`\tau _{\mathrm{sca}}`$ of the two-phase models models are similar to those of the hierarchical models with uniform dust. The contrast of both with purely hierarchical models illustrates the importance of dust-free regions (in real space, possibly caused by extensions of hot, low-density material into the nebulae).

Figure 4 shows what we meant by saying that Figure 1 was produced by a typical hierarchical model – one with a typical $`\tau _{\mathrm{ext}}`$. Those with large $`\tau _{\mathrm{ext}}`$ have the star embedded within dusty material, and a relatively low central dust density leads to a low $`\tau _{\mathrm{ext}}`$.

Differences among models, and errors arising from interpreting reflection nebulae with uniform models, increase with optical depth. Figure 5 shows the same as Figure 4 ($`a`$ = 0.6; $`g`$ = 0.6), except that $`\tau _0`$= 4 instead of 2. The open squares are for various initial dust distributions with $`D`$ =2.3, with 33% of the density uniform; the filled squares, the same with $`D`$ = 2.6; open triangles, purely hierarchical density, $`D`$=2.3; filled triangles, purely hierarchical with $`D`$ = 2.6. The circles enclose the entire range of our two-phase models models.

Figure 5 shows an important result: optically thick hierarchical models can show very modest values of $`\tau _{\mathrm{ext}}`$, but interpreting the nebula with uniform models will likely result in a severe underestimate of the albedo unless there is material between the clumps to scatter light efficiently. For the purely hierarchical models with $`D`$ = 2.3 (open triangles), there is very little scattered light between the clumps (i.e., large $`\tau _{\mathrm{sca}}`$) and limited stellar extinction, resulting in albedos predicted from uniform models that are as low as 0.17. The spread among the two-phase models (circles) is relatively small because at large optical depth they become statistically uniform. Since they are clumpy, the uniform-model albedos they suggest are significantly lower than the 0.6 used in the actual two-phase models.

The squares in Figure 5 have 33% of the dust mass in a constant density, so that the minimum optical depth is 4/3. This addition greatly increases the scattered flux and makes the difference between $`D`$ = 2.3 and 2.6 almost unimportant. The albedo that would be derived from uniform models is increased to $``$0.5 over the purely hierarchical models, still short of the correct value of 0.6.

The spread of the $`\tau _{\mathrm{ext}}`$ and $`\tau _{\mathrm{sca}}`$ in Figures 4 and 5 show that the specification of $`\tau _0`$ and $`D`$ (along with grain properties) does not specify the scattered and stellar fluxes from reflection nebulae in even a statistical sense. By contrast, $`D`$ constrains the statistical properties of the projected brightness distribution for both line and continuum emission from the gas. The large variations in properties of reflection nebulae with viewing angle come about because the radiation arises from the geometry of the statistical distribution of the ISM relative to a point source rather than to itself. Emissions form gas involve an integration over larger regions of more equal weight.

It is unrealistic to suppose that reflection nebulae are spherical with a central exciting star. We have considered only this geometry because we feel that it provides rather general results for strongly forward-throwing scattering (say, $`g0.5`$). In this case, the dust on the near side of the nebula is likely to provide most of the scattered light. We consider our models as representing only this near-side dust.

## 4 Discussion

### 4.1 Reflection Nebulae

Our hierarchical models, with their overall spherical distribution of clumps, are primarily aimed at interpreting reflected light from reflection nebulae illuminated by an exciting star. We will discuss reflection of the diffuse Galactic radiation field from dusty globules at the end of this section.

We now examine the limits on $`a`$ imposed by the amount of scattered light observed in reflection nebulae. Since real nebulae show the illuminated cloud as it is viewed from only one direction, we must consider the spread of ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) from any particular dust distribution. We seek $`a_{\mathrm{max}}`$ and $`a_{\mathrm{min}}`$, the minimum and maximum values of $`a`$ that can fit the observations.

For illustrative purposes, we consider the 1300 Å observations of NGC 7023 with 50% of the stellar extinction assumed to be within the nebula (Witt et al. 1993; the circled dot in Figures 1 and 2). In Figure 1 it lies in the middle of the envelope of the ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) points from various viewing angles. If we decreased the $`a`$ of the model, the scattered flux at each direction would become fainter and $`\tau _{\mathrm{sca}}`$ larger. The stellar fluxes would be unchanged because they depend only on extinction, not scattering. Thus, the pattern of points would move up in the diagram, and the point could no longer be fitted. To find the $`a_{\mathrm{min}}`$ for NGC 7023, we must vary $`\tau _0`$ and initial seeds (dust distributions) and find the minimum $`a`$ for which the envelope of individual points lies just above the observation. A similar procedure leads to $`a_{\mathrm{max}}`$. Of course, it is improbable that the actual $`a`$ is as low as $`a_{\mathrm{min}}`$ because the observed ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) would probably not lie at the edge of the ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) distribution as seen from a random direction.

We have performed this exercise and found the following:

(a) All of our models required an albedo of $`0.5`$ to fit the lower point of NGC 7023. Varying parameters and geometries to obtain the absolute minimum $`a`$ (0.46; see below) produced only changes in details. This result is robust because the observed $`\tau _{\mathrm{ext}}`$ is near the most efficient value for producing scattered light, as suggested by the uniform models (see the minima of $`\tau _{\mathrm{sca}}`$ in the curves in Figure 1).

(b) The scattered light is near its maximum when $`\tau _0`$ $``$2, where the mean $`\tau _{\mathrm{ext}}`$ is about unity, so $`a_{\mathrm{min}}`$ is best estimated for our hierarchical models near this value. The maximum scattered light for uniform models occurs at $`\tau _0`$ $``$ 1.1 – 1.4, increasing slowly with the albedo.

(c) The smallest $`a_{\mathrm{min}}`$ we found for the circled dot in Figure 1 is 0.48. It occurs for $`\tau _0`$ $``$2, $`D`$ = 2.6, and a particular seed whose ($`\tau _{\mathrm{ext}}`$, $`\tau _{\mathrm{sca}}`$) point is low in Figure 4. The effects of $`g`$ are minimal on the locus of the lower envelope of the individual points; the minimum $`a`$ is achieved for $`g0.5`$. We tried such a wide range of models that we doubt that there is an $`a_{\mathrm{min}}`$ drastically lower.

(d) No useful estimate of $`a_{\mathrm{max}}`$ is possible because inefficient scattering can be produced by extreme clumping (as small a $`D`$ as allowed) and with no uniform component. For $`D`$ = 2.3, an $`a1`$ is needed to fit NGC 7023 for some distributions of clumps. For $`D`$ = 2.6, $`a_{\mathrm{max}}`$ is about 0.8. The dust geometry that produced $`a_{\mathrm{min}}`$ did not allow $`a`$ to be above 0.6, but it was chosen to have especially efficient scattering. Alternatively, very optically thick models that have low stellar extinctions as seen from some viewing angles can fit observed fluxes with very high albedos. For instance, we have fitted the NGC 7023 point in Figure 1 with two-phase models, $`\tau _0`$ = 8, $`a`$ = 0.85. We have no doubt that many optically thick hierarchical models could fit the observations as well. Perhaps such extreme models could be ruled out with far-infrared (FIR) fluxes, although the low absorption tends to compensate for the high optical depth. Furthermore, the FIR arises from dust surrounding the star in all directions, while the scattering is mainly from dust close to the line of sight. In any case, these models are hardly reasonable on physical grounds.

Voyager 2 and HUT observed NGC 7023 at 1000Å (Witt et al. 1993). Uniform models predict $`a`$0.45 if 40% of the stellar extinction occurs within the nebula. Hierarchical models can fit the observation (not plotted, but in the upper right corner of Figures 1 – 3) with $`a`$ = 0.8 with $`\tau _04`$. Of course, we make no claim that such a high $`a`$ is correct. On the other hand, we can fit the point with $`a`$ = 0.4, lower than the uniform models. Once again we see that the observations are subject to a very ambiguous interpretation.

Burgh, McCandliss, & Feldman (2002) used a uniform model for NGC 2023 with $`\tau _{\mathrm{ext}}`$ = 2.4. They derived an albedo of $`0.39_{0.05}^{+0.12}`$ at 1345 Å. Our hierarchical model shown in Figure 1, with $`a`$ = 0.6, could accommodate their observations, and others that scatter light inefficiently could have $`a1`$.

Figure 6 shows the albedos derived from uniform models, using the scattering and extinctions from the hierarchical model $`(a,g,D,\tau _0)`$ = (0.6, 0.6, 2.3, 1), plotted against the stellar extinctions at each viewing angle. The albedo used to generate the scattered fluxes and extinctions, 0.6, is marked by the dashed line. The distribution of points for the hierarchical model with $`\tau _0`$ = 2 is virtually identical. We see that the minimum $`a`$ derived from uniform models is $``$0.32; the average for $`0.6\tau _{\mathrm{ext}}1.6`$ is $``$0.42. At moderate extinction, leading to bright scattered flux, uniform models underestimate $`a`$. At small $`\tau _{\mathrm{ext}}`$ uniform models overestimate $`a`$. A hierarchical model in which the star happens to be within a clump has smaller errors. Adding a constant density decreases the errors at $`\tau _{\mathrm{ext}}`$ $``$ 1 so that the average $`a`$ is $``$ 0.5. Of course, there is no error in the unrealistic case of the density being uniform throughout.

An important question is, how reliably can we determine $`a(\lambda )`$ from observations of a particular nebula covering a significant range in wavelengths, implying differing opacities or optical depths? In this case, we might hope to find relative optical properties of the grains accurately even with uniform models. Figure 7 shows the results for two wavelengths for which the extinctions (both $`\tau _0`$ and $`\tau _{\mathrm{ext}}`$) differ by a factor of 2: for example, (1000 Å)/(1540 Å) if the extinction law is the Galactic average (Cardelli, Clayton, & Mathis 1989) and $`R_V=A(V)/[A(B)A(V)]=3.1`$. The figure displays the ratio of $`a`$($`\tau _0`$=2) to $`a`$($`\tau _0`$=1), both albedos derived from uniform models, for the same hierarchical model as shown in Figure 1. The abscissa in the figure is $`\tau _{\mathrm{ext}}`$ for the lower extinction. All points in the figure were computed with the $`a`$ = 0.6, so the “true” value of the albedo ratio is unity. We see that uniform models have difficulty in predicting even relative values of $`a(\lambda )`$, since the spread in $`a(\tau _0=2)/a(\tau _0=1)`$ is $``$40%. The figure gives the impression that if $`\tau _{\mathrm{ext}}`$ is $`1.5`$, uniform models would derive a greater albedo at short wavelengths (say, 1000 Å) than at longer (1540 Å for a factor of 2 in the extinction). For NGC 7023 the opposite is found to be the case (Murthy et al. 1993; Witt et al. 1993). From Figure 7 we would conclude that $`a`$(1000 Å)$`<a`$(1540 Å) is robust, since Witt et al. (1993) suggest $`\tau _{\mathrm{ext}}`$(1540 Å) $``$(1.5 – 2). However, models with the star embedded in a clump have the albedo ratios shift from mostly $`<1`$ to mostly $`>1`$ at $`\tau _{\mathrm{ext}}`$(1540 Å) $`2.5`$, considerably greater than shown in the figure. These models would produce an illusory result $`a`$(1000 Å) $`<a`$(1740 Å).

The basic problem is that the importance of a clump depends upon its optical depth. If a clump becomes optically thin, its existence no longer very significant to the radiative transport. Similarly, increasing the optical depth of an optically thick clump has a limited effect on the radiation. While the density structure within the ISM is independent of the wavelength, its geometry as regards radiation transfer is not. Since the clumping produces the large spread in $`\tau _{\mathrm{sca}}`$($`\theta ,\varphi )`$, it is not surprising that there is a wide spread in the importance of the clumps among various viewing angles.

The problems of determining the variation of albedo at various wavelengths are small if the ratio of optical depths is close to unity. Statements that the albedo is smaller at the 2175 Å feature than on either side of it (e.g., Witt et al. 1992) seem robust, besides being completely plausible physically. At nearby wavelengths $`A(\lambda )/N(`$H), the extinction per H nucleus, can be the same at three wavelengths, and comparison of the relative albedos between those wavelengths is completely safe.

We have seen that several properties of reflection nebulae change if the exciting star is embedded within a rather dense clump of dust. There is an increase in the scattered light; the lowest values of $`\tau _{\mathrm{sca}}`$ for purely hierarchical models in Figure 4$`a`$ (i.e., the large amounts of scattered light for the purely hierarchical models) belong to such cases. The central clump ensures that the angle-averaged extinction, $`\tau _{\mathrm{ext}}`$, is larger than for most models with the same $`\tau _0`$ (i.e., the points in Figure 4$`a`$ tend to have large $`\tau _{\mathrm{ext}}`$ if $`\tau _{\mathrm{sca}}`$ is small). Figure 4$`b`$ shows that the addition of a constant density component makes the effects of a central clump rather unimportant because the radiation scattered by the uniform dust dominates. Another effect of a star being embedded in a dense clump is that the errors of the albedo derived from observations of the same object at different wavelengths, illustrated in Figure 6, are reduced because more of the scattering arises from the central clump. For such a clump, the opacity, but not the geometry, changes with wavelength. Finally, there would tend to be more FIR radiation from an embedded star because the grains are relatively close to the star, so they are absorbing a relatively large stellar flux, and the mean temperature of the grains would be relatively large for the same reason.

Galactic radiation scattered by isolated globules seen away from the Galactic plane provides another diagnostic of grains (Mattila 1980). Witt, Oliveri, & Schild (1990) have analyzed such a globule and give references to other examples and discussion of the process. The density is concentrated towards the center of such globules. The reflected intensity rises outward from the center, peaks, and falls off to the sky brightness at the edge. showing that the scattering is forward-throwing. The reflected intensity can be related to the incident radiation predicted for both stars and the DGL out of the Galactic plane. Results were that at 4700 Å grains are forward-throwing ($`g`$0.8) and very efficient at scattering ($`a0.8`$). The globule was modeled with a very centrally condensed density distribution that is smooth and spherically symmetric, as seemed completely appropriate at the time. Since turbulence driving clumps in the diffuse ISM may decay when driving into a steepening density gradient, perhaps the ISM within globules is relatively free from clumping. We cannot assess the uncertainties in the grain properties derived from globules until we know more about the clumping within them. The major uncertainty would then be the prediction of the radiation incident upon them from the hierarchically clumped Galaxy.

### 4.2 The Diffuse Galactic Light

There is no doubt about the applicability of hierarchical clumping to the DGL. Observations such as those at 1740 Å (Schiminovich et al. 2001) are based on large portions of the sky. We have seen that uniform spherical models of reflection nebulae underestimate $`a`$ if applied to optically thick situations (large $`\tau _0`$, but not necessarily large $`\tau _{\mathrm{ext}}`$). Our reflection nebula models are not ideally suited for interpreting the DGL in the UV, since distributions of the dust and exciting stars are skewed towards the Galactic plane. However, our results have strong implications as regards the interpretation of the DGL as discussed in Henry (2002, hereafter H02). In that paper it was concluded that $`a0.1`$ for $`\lambda 1400`$ Å.

H02 performed radiative transport calculations by using $`B,V`$ magnitudes and the spectral type for each star in the entire Hipparchos Input Catalog. The reddening, $`E(BV)`$, follows from $`(BV)`$ and spectral type. The extinction $`A_\lambda `$ is $`E(BV)R_V(A_\lambda /A_V)`$, where $`R_V`$ is usually taken to be 3.1. The $`\tau _{\mathrm{ext}}`$ for each star is $`A(\lambda )`$/1.086.

The heart of H02 is its unique treatment of radiative transfer. The scattering is taken to be

$$\mathrm{exp}(\tau _{\mathrm{sca}})=a(1e^{\tau _{\mathrm{ext}}})\mathrm{exp}\{(1a)[\tau _{\mathrm{ext}}+\mathrm{ln}(1+e^{\tau _{\mathrm{ext}}})\mathrm{ln}2]\},$$
(2)
independently of the geometry. However, Witt, Friedmann, & Sasseen (1997) have shown from clumpy models that the distribution of cloud optical depths can strongly affect the interstellar radiation field.

At low optical depths, the H02 relation between $`\tau _{\mathrm{ext}}`$ and $`\tau _{\mathrm{sca}}`$ approaches the uniform model. The dotted lines in Figure 1 show the H02 models for three albedos identified by the nearby solid lines. We see that the H02 assumption predicts more scattered light than uniform models, and a lower albedo follows.

For $`\lambda `$$`2000`$ Å, the $`\tau _{\mathrm{ext}}`$ for most hot stars is probably rather large, since the mean $`A(V)`$ is $``$1 magnitude kpc-1 and the ultraviolet is 3 – 5 times larger. For $`\tau _{\mathrm{ext}}`$ = 1.5, Figure 1 shows that the H02 model with $`a`$ = 0.6 predicts a scattered flux $``$exp(0.4) = 1.5 times larger than that “typical” hierarchical model, chosen for illustrative purposes only. The conclusion of H02 that $`a`$(1000 Å) is very low (“perhaps 0.1”) is perhaps understandable. Such a low albedo is not compatible with the observations of NGC 7023 by Witt et al. (1993), in which the scattered flux was larger than the stellar, or Murthy et al. (1993), which suggests that $`a`$(1150 Å) $`0.55`$ if 40% of the extinction is within the nebula, as assumed by Witt et al. (1992).

Murthy, Henry, & Holberg (1991) used Voyager 2 to search for DGL for $`\lambda 1300`$ Å. They failed to detect any, possibly confirming the H02 picture. Our reflection nebula results commonly show that the same nebula can show very low and rather high scattered fluxes as seen from different angles. We strongly suspect that the same property applies to the DGL. Our general result from hierarchical models is that photons, especially the unscattered ones coming directly from the star, can leak out of dust distributions far more easily than either uniform-density or H02 models would predict, so a dark sky does not necessarily imply a low albedo.

## 5 Final Remarks and Summary

We have arrived at a pessimistic, but we feel realistic, assessment of the ability of reflection nebulae illuminated by embedded stars to serve as diagnostics of interstellar grain properties. It is better to realize what we do not know, rather than believe that we know more than we do! Even if the scattered light from a reflection nebula and the extinction suffered by its exciting star are observed perfectly, the geometry of the dust is vital in the interpretation of the observations. Hierarchical clumping, “clumps within clumps within clumps”, has a very large impact on the interpretation of both reflection nebulae and the DGL as compared to smooth density structures. Hierarchical structures, demonstrated by CO, 100 $`\mu `$m, various other molecules, and H I observations, make any determination of the optical properties of grains from reflection nebulae very difficult. Our models had plausible power spectra for the autocorrelation function of the projected densities, with an exponent of the power of $`3`$, unlike either uniform models or those that are clumped on a single scale (exponent $``$0). The exact geometry of the dust is very difficult to determine because it depends on the placement of hierarchical clumps of dust around the star.

We have made no interpretations of any reflection nebulae. It is now relatively easy to use Monte Carlo transfer codes to predict the scattering (and polarization) for objects with a known geometry. Inverting observed intensities to determine the correct geometry, though, is very difficult, even with observations with good spatial resolution taken in multiple wavelengths. Our hierarchical models have been schematic and do not begin to exhaust the possibilities of the real geometry of any reflection nebula. Observations in the infrared and line emissions with good spatial resolution may eventually put useful limits on the actual geometry of the gas and dust, but the problem of inverting the data will always be difficult.

An interesting possibility is that isolated Bok globules may not be strongly clumped, so reflection of Galactic light from them can diagnose grain properties. The penetration of turbulence into them might be weaker than its percolation through larger structures. Nonclumpy globules can serve to diagnose grains because they scatter incident Galactic radiation (see Witt et al. 1990 for references and an explanation).

Each of our hierarchical models shows a very wide range of extinction optical depths ($`\tau _{\mathrm{ext}}`$) and scattering ($`\tau _{\mathrm{sca}}`$) when the nebula is viewed from various directions. The $`\tau _{\mathrm{ext}}`$ and $`\tau _{\mathrm{sca}}`$ observed for a reflection nebula represent a view of the dust from one direction. There is a very wide range of optical properties ($`a,g`$) that can fit a given observation. With one or more hierarchical models we can fit the best observed UV reflection nebula, NGC 7023, with albedos $``$ 0.5 over the range of $`g`$ that we tested (0.6 – 0.85) and with a variety of averaged optical depths through the model.

In general, reflection nebulae have substantial optical depths ($`\tau _{\mathrm{ext}}`$ $`13`$); otherwise, the scattered fluxes are too faint. At substantial optical depths, uniform models can underestimate the albedo rather severely because they overestimate the scattered light. Unless the star happens to be embedded within a dense clump, radiation can leak out of hierarchical clumps much more easily than from uniform dust.

The Diffuse Galactic Light (DGL) has been interpreted (Henry 2002) with a recipe that predicts much more scattering from each star than our models, resulting in a low estimates of the UV albedo. Murthy et al. (1991) failed to detect UV DGL, but our fractal models often have many directions from which the star and scattered light are very faint. Thus, we do not believe that the faintness of the DGL in particular directions necessarily signifies a low albedo.

The determination of both $`a`$ and $`g`$ depends upon the variation of reflected intensity across the face of a centrally illuminated reflection nebula, while we have only considered the flux of the scattered light. We feel that if the flux is as weakly constrained as our models show, the intensity will be similarly subject to variation because of the unknown placement of the star relative to the surrounding dust.

Even the relative variation of albedo with wavelength is difficult to determine from reflection nebulae if there are significant changes in opacity among the wavelengths. The problem is that the geometry is not really the same when the opacity changes. As opacity per H atom increases, optically thin clumps become thick and scatter light with a different geometrical arrangement, so the geometry of the nebula depends on wavelength.

Our overall assessment is that the optical properties of grains are probably as well constrained by theory as by observations. This statement is made in spite of well-known uncertainties in the theory of interstellar grains. These uncertainties are major: whether typical large grains are chemically homogeneous (e.g., silicate or carbonaceous) or composite, whether grains contain voids, or if grains have very loose (“fractal”) structure. However, theory does not permit the large range of possible albedos provided by reflection nebulae. At first glance, various theories seem quite different. Weingartner & Draine (2001) have PAH molecules, small grains, and large grains of silicate or graphite. Li & Greenberg (1997) have PAHs, small grains, and large grains with silicate cores and organic refractory mantles. Mathis (1996) has composite grains with silicates, carbon, and vacuum intermixed within the same grain, plus small graphite that mimics PAHs. Draine & Lee (1984) have silicate and graphite with a truncated power-law distribution of sizes. Each size distribution is different. They all predict roughly the same $`a(\lambda )`$ because the response of the grains to radiation is set by fitting the well-observed interstellar extinction law $`\tau _{\mathrm{ext}}`$($`\lambda `$), and cosmic abundances (though somewhat controversial themselves) constrain the materials. For instance, the values predicted for $`a`$(1430 Å) are (0.40, 0.30, 0.37, 0.40) for Weingartner & Draine (2001), Li & Greenberg (1997), Mathis (1996), and Draine & Lee (1984), respectively. The $`A(V)`$ are (0.65, 0.63, 0.53, 0.57).

Fortunately, optical depths and column densities derived from comparisons of absorption and scattering along the same very narrow beam (i.e., towards a star) are not affected by the microstructure within the ISM. These measurements include most determinations of depletions of various ions in the ISM, as well as the wavelength dependence of the interstellar extinction law. Other implications for the hierarchical nature of dust clumping include a greatly increased penetration of radiation into molecular clouds and PhotoDissociation Regions (PDRs).

We would like to thank Prof. A. Witt for constructive comments on a draft of this paper that improved it considerably. BW and KW acknowledge financial support from NASA’s Long Term Space Astrophysics Research Program, NAG5 6039 (KW), NAG5 8412 (BW), the National Science Foundation, AST 9909966 (BW and KW), and a PPARC Advanced Fellowship (KW). We appreciate the constructive and prompt comments of two anonymous referees. This research has benefited from the NASA Astrophysics Data System and the SIMBAD database.


text:
# Spinor Bose-Einstein Condensates with Many Vortices

## Abstract

Vortex-lattice structures of antiferromagnetic spinor Bose-Einstein condensates with hyperfine spin $`F=1`$ are investigated theoretically based on the Ginzburg-Pitaevskii equations near $`T_c`$. The Abrikosov lattice with clear core regions are found never stable at any rotation drive $`\mathrm{\Omega }`$. Instead, each component $`\mathrm{\Psi }_i`$ $`(i=0,\pm 1)`$ prefers to shift the core locations from the others to realize almost uniform order-parameter amplitude with complicated magnetic-moment configurations. This system is characterized by many competing metastable structures so that quite a variety of vortices may be realized with a small change in external parameters.

Realizations of the Bose-Einstein condensation (BEC) in atomic gases have opened up a novel research field in quantized vortices as created recently with various techniques JILA99 ; ENS00 ; MIT01 ; JILA01 . Especially interesting in these systems are vortices of spinor BEC’s in optically trapped 23Na Ketterle98  and 87Rb Barrett01 , where new physics absent in superconductors Tinkham , 4He Donnelly , and 3He Volovik87 ; Thuneberg99 ; Kita01 , may be found.

Theoretical investigations on spinor BEC’s were started by Ohmi and Machida OM98  and Ho Ho98 , followed by detailed studies on vortices with a single circulation quantum Yip99 ; Busch99 ; Volovik00 ; Marzlin00 ; Khawaja01 ; Martikainen01 ; Kurihara01 ; Isoshima01 ; Isoshima02 ; Mizushima02 . However, no calculations have been performed yet on structures in rapid rotation where the trap potential will play a less important role. Indeed, the clear hexagonal-lattice image of magnetically trapped 23Na MIT01  suggests that predictions on infinite systems are more appropriate for BEC’s with many vortices. Such calculations have been carried out by Ho for the single-component BEC Ho01  and by Mueller and Ho for a two-component BEC Ho02  near the upper critical angular velocity $`\mathrm{\Omega }_{c2}`$ at $`T=0`$.

The purpose of the present paper is to perform detailed calculations on vortices of $`F=1`$ spinor BEC’s in rapid rotation to clarify their essential features. To this end, we focus on the mean-field high-density phase rather than the low-density correlated liquid phase Cooper01 , and use the phenomenological Ginzburg-Pitaevskii (or Ginzburg-Landau) equations near $`T_c`$ GP58 ; GA82  instead of the Gross-Pitaevskii equations at $`T=0`$. Since fluctuations are small in the present system, this approach will yield quantitatively correct results near $`T_c`$. It should be noted that the corresponding free-energy is formally equivalent to that derived with Ho’s “mean-field quantum Hall regime” near $`\mathrm{\Omega }_{c2}`$ at $`T=0`$ Ho01 ; Ho02 , so that the results obtained here are also applicable to that region.

Model.— The free-energy density of an $`F=1`$ spinor BEC near $`T_c`$ may be expanded with respect to the order parameters $`\mathrm{\Psi }_i`$ $`(i=0,\pm 1)`$ as

$`f=\alpha \mathrm{\Psi }_i^{}\mathrm{\Psi }_i+\mathrm{\Psi }_i^{}{\displaystyle \frac{\left(i\mathrm{}\mathbf{}M𝛀\times 𝐫\right)^2}{2M}}\mathrm{\Psi }_i`$
$`+{\displaystyle \frac{\beta _\mathrm{n}}{2}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}\mathrm{\Psi }_i\mathrm{\Psi }_j+{\displaystyle \frac{\beta _\mathrm{s}}{2}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}(F_\mu )_{ik}(F_\mu )_{jl}\mathrm{\Psi }_k\mathrm{\Psi }_l.`$ (1)
Here $`\alpha `$, $`\beta _\mathrm{n}`$ , and $`\beta _\mathrm{s}`$ are expansion parameters, $`F_\mu `$ $`(\mu =x,y,z)`$ denotes the spin operator, $`M`$ is the particle mass, and summations over repeated indices are implied. The rotation axis $`𝛀`$ is taken along $`𝐳`$. The quantities $`\beta _\mathrm{n}`$ and $`\beta _\mathrm{s}`$ are assumed to be constant near $`T_c`$ with $`\beta _\mathrm{n}>0`$, whereas $`\alpha `$ changes its sign at $`T_c`$ with $`\alpha >0`$ for $`T<T_c`$. To simplify Eq. (1), we measure the length, the energy density, the angular velocity, and the order parameter in units of $`\mathrm{}/\sqrt{2M\alpha }`$, $`\alpha ^2/\beta _\mathrm{n}`$, $`\alpha /\mathrm{}`$, and $`\alpha /\beta _\mathrm{n}`$, respectively. The corresponding free-energy density is obtained from Eq. (1) by $`\alpha 1`$, $`\beta _\mathrm{n}1`$, $`\mathrm{}^2/2M1`$, $`M𝛀\frac{1}{2}𝛀`$, and $`\beta _\mathrm{s}g_\mathrm{s}\beta _\mathrm{s}/\beta _\mathrm{n}`$. It thus takes a simple form with only two parameters $`(g_\mathrm{s},\mathrm{\Omega })`$. We then introduce a couple of operators by $`a\mathrm{}(_x+i_y)/\sqrt{2}`$ and $`a^{}\mathrm{}(_x+i_y)/\sqrt{2}`$ with $`\mathbf{}\mathbf{}\frac{i}{2}𝛀\times 𝐫`$ and $`\mathrm{}1/\sqrt{\mathrm{\Omega }}`$ which satisfy $`aa^{}a^{}a=1`$. Equation (1) now reads

$`f=`$ $`\mathrm{\Psi }_i^{}[(2a^{}a+1)\mathrm{\Omega }1]\mathrm{\Psi }_i+{\displaystyle \frac{1}{2}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}\mathrm{\Psi }_i\mathrm{\Psi }_j`$ (2)
$`+{\displaystyle \frac{g_\mathrm{s}}{2}}\mathrm{\Psi }_i^{}\mathrm{\Psi }_j^{}(F_\mu )_{ik}(F_\mu )_{jl}\mathrm{\Psi }_k\mathrm{\Psi }_l,`$
with $`\mathrm{\Omega }_{c2}=1`$, and the free-energy is given by

$`{\displaystyle f(𝐫)𝑑𝐫}.`$ (3)
We can find the stable structure for each $`(g_\mathrm{s},\mathrm{\Omega })`$ by minimizing $``$. We have performed extensive calculations over the whole antiferromagnetic region $`g_\mathrm{s}0`$, where $`(\mathrm{\Psi }_1,\mathrm{\Psi }_0,\mathrm{\Psi }_1)=\mathrm{e}^{i\theta }𝒰(0,1,0)`$ and $`f=\frac{1}{2}`$ at $`\mathrm{\Omega }=0`$ with $`\theta `$ an arbitrary phase and $`𝒰`$ the spin-space rotation Ho98 .

A major difference from superfluid 3He Kita01  lies in the fact that terms such as $`\mathrm{\Psi }_i^{}aa\mathrm{\Psi }_j`$ $`(ij)`$ are absent, i.e., there are no gradient couplings between different components. Hence $`\mathrm{\Omega }_{c2}`$ is the same for all components, whereas in 3He only a single component becomes finite at $`\mathrm{\Omega }_{c2}`$ to realize the polar state Kita01 . This degenerate feature is what characterizes the present system to bring many competing metastable structures, as seen below.

Method.— We minimize Eq. (3) with the Landau-level-expansion method (LLX) Kita98 ; Kita01  by expanding the order parameters as

$$\mathrm{\Psi }_i(𝐫)=\sqrt{𝒱}\underset{N=0}{\overset{\mathrm{}}{}}\underset{𝐪}{}c_{N𝐪}^{(i)}\psi _{N𝐪}(𝐫),$$
(4)
with $`N`$ the Landau-level index, $`𝐪`$ the magnetic Bloch vector, and $`𝒱`$ the system volume. The basis functions $`\{\psi _{N𝐪}\}`$ are eigenstates of the magnetic translation operator $`T_𝐑\mathrm{exp}[𝐑(\mathbf{}+\frac{i}{2}𝛀\times 𝐫)]`$, which can describe periodic vortex structures of all kinds Kita98 . Its explicit expression is given by

$`\psi _{N𝐪}(𝐫)=`$ $`{\displaystyle \underset{n=𝒩_\mathrm{f}/2+1}{\overset{𝒩_\mathrm{f}/2}{}}}\mathrm{e}^{i[q_y(y+\mathrm{}^2q_x/2)+na_{1x}(y+\mathrm{}^2q_xna_{1y}/2)/\mathrm{}^2]}`$ (5)
$`\times \mathrm{e}^{ixy/2\mathrm{}^2(x\mathrm{}^2q_yna_{1x})^2/2\mathrm{}^2}`$
$`\times \sqrt{{\displaystyle \frac{2\pi \mathrm{}/a_2}{2^NN!\sqrt{\pi }𝒱}}}H_N({\displaystyle \frac{x\mathrm{}^2q_yna_{1x}}{\mathrm{}}}),`$
with $`𝒩_\mathrm{f}^2`$ the number of the circulation quantum $`\kappa h/M`$ in the system, $`H_N`$ the Hermite polynomial, and $`𝐚_j`$ the basic vectors in the $`xy`$ plane with $`𝐚_2\widehat{𝐲}`$ and $`a_{1x}a_2=2\pi \mathrm{}^2`$. Substituting Eq. (4) into Eq. (3) and carrying out the integration in terms of $`(s,t)(x/a_{1x},y/a_2xa_{1y}/a_{1x}a_2)`$, the free energy is transformed into a functional of the expansion coefficients $`\{c_{N𝐪}^{(i)}\}`$, the apex angle $`\beta \mathrm{cos}^1(a_{1y}/a_1)`$, and the ratio of the two basic vectors $`\rho a_2/a_1`$ as $`=[\{c_{N𝐪}^{(i)}\},\beta ,\rho ]`$. For a given $`\mathrm{\Omega }`$, we directly minimize $`/𝒱`$ with respect to these quantities.

Search for stable structures.— We here sketch our strategy to find stable structures. To this end, we first summarize the basic features of the conventional Abrikosov lattice within the framework of LLX Kita98 : (i) any single $`𝐪=𝐪_1`$ suffices, due to the broken translational symmetry of the vortex lattice; (ii) the triangular (square) lattice is made up of $`N=6n`$ ($`4n`$) Landau levels ($`n=0,1,2\mathrm{}`$); (iii) more general structures can be described by $`N=2n`$ levels, odd $`N`$’s never mixing up since those bases have finite amplitudes at the cores. This Abrikosov lattice has a single circulation quantum per unit cell.

With multi-component order-parameters, there can be a wide variety of vortices, which may be divided into two categories. We call the first category as “shift-core” states, where core locations are different among the three components with an enlarged unit cell. General structures with $`n_\kappa `$ circulation quanta per unit cell can be described by using $`n_\kappa `$ different $`𝐪`$’s, where the unit cell becomes $`n_\kappa `$ times as large as that of the Abrikosov lattice. For example, structures with two quanta per unit cell are given by choosing $`(𝐪_1,𝐪_2)=(\mathrm{𝟎},\frac{𝐛_1+𝐛_2}{2})`$, where $`𝐛_1`$ and $`𝐛_2`$ are reciprocal lattice vectors. This possibility was already considered by Mueller and Ho Ho02  for a two-component system and shown to yield stable structures. It also describes the mixed-twist lattice to be found in 3He Kita01 . The second category may be called “fill-core” states with a single circulation quantum per unit cell (i.e., a single $`𝐪`$ is relevant). Here the cores of the conventional Abrikosov lattice are filled in by some superfluid components using odd-$`N`$ wavefunctions of Eq. (5). This entry of odd-$`N`$ Landau levels occurs as a second-order transition below some critical angular velocity smaller than $`\mathrm{\Omega }_{c2}/3`$. It has been shown that the $`A`$-phase-core vortex in the $`B`$-phase of superfluid 3He belong to this category Kita01 .

We have carried out an extensive search for stable structures with up to $`n_\kappa =9`$ circulation quanta per unit cell, including fill-core states. Since we are near $`T_c`$ where normal particles are also present, we have performed the minimization without specifying the value of the magnetic moment $`M`$ for the superfluid components. However, all the stable states found below have $`M=0`$. Each of the three components were expanded as Eq. (4) using $`n_\kappa `$ different $`𝐪`$’s, and the free-energy (3) is minimized with respect to $`c_{N𝐪}^{(i)}`$, $`\beta `$, and $`\rho `$ by using Powell’s method NR . To pick out stable structures correctly, we calculated Eq. (3) many times starting from different initial values for $`c_{0𝐪}^{(i)}`$, $`\beta `$, and $`\rho `$ given randomly within $`0.1\mathrm{Re}c_{0𝐪}^{(i)},\mathrm{Im}c_{0𝐪}^{(i)}0.9`$, $`0.1\pi \beta 0.5\pi `$, and $`0.8\rho 3`$, respectively. The state with the lowest energy was thereby identified as the stable structure. The spin quantization axis and an overall phase were fixed conveniently to perform efficient calculations. Thus, any structures obtained from the solutions below with a spin-space rotation and a gauge transformation are also stable.

Instability of Abrikosov lattice.— The present calculations have revealed that the Abrikosov lattice with clear core regions is never stable at any rotation drive $`\mathrm{\Omega }`$ over the entire antiferromagnetic region $`g_\mathrm{s}0`$. Thus, any optical experiments to detect vortices by amplitude reductions are not suitable for the spinor vortices.

The most stable Abrikosov lattice is given by

$`\mathrm{\Psi }_0(𝐫)=\sqrt{V}{\displaystyle \underset{N}{}}c_N^{(0)}\psi _{N𝐪_1}(𝐫),`$ (6)
with $`c_{N𝐪_1}^{(0)}`$ real, $`N=6n`$, $`\beta =\frac{\pi }{3}`$, and $`\rho =1`$. Here the antiferromagnetic component $`\mathrm{\Psi }_0`$ forms a hexagonal lattice with a single circulation quantum per unit cell. Below some critical velocity $`\mathrm{\Omega }_\mathrm{f}`$ smaller than $`\mathrm{\Omega }_{c2}/3`$, the core regions start to be filled in by

$`\mathrm{\Psi }_{\pm 1}(𝐫)=i\sqrt{V}{\displaystyle \underset{N^{}}{}}c_N^{}^{(1)}\psi _{N^{}𝐪_1}(𝐫),`$ (7)
with $`N^{}`$ odd. The second transition for this odd-Landau-level entry occurs at $`\mathrm{\Omega }_\mathrm{f}=0.1497\mathrm{\Omega }_{c2}`$ and $`0.0938\mathrm{\Omega }_{c2}`$ for $`g_\mathrm{s}=0.1`$ and $`0.3`$, respectively.

However, calculations down to $`0.0001\mathrm{\Omega }_{c2}`$ of using $`1800`$ Landau levels for $`g_\mathrm{s}0`$ have clarified that the above fill-core state carries higher free energy than the following shift-core state with two circulation quanta per unit cell:

$`\mathrm{\Psi }_1(𝐫)=\sqrt{V}{\displaystyle \underset{N}{}}c_N\psi _{N𝐪_1}(𝐫),`$ (8)
$`\mathrm{\Psi }_1(𝐫)=\sqrt{V}{\displaystyle \underset{N}{}}c_N\psi _{N𝐪_2}(𝐫),`$ (9)
with $`c_N`$ real and common to both, $`N=4n`$, $`\beta =\frac{\pi }{2}`$, $`\rho =1`$, and $`(𝐪_1,𝐪_2)=(\mathrm{𝟎},\frac{𝐛_1+𝐛_2}{2})`$. The cores of $`\mathrm{\Psi }_{\pm 1}(𝐫)`$ are shifted from each other by $`\frac{1}{2}(𝐚_1+𝐚_2)`$. Figure 1 displays basic features of this shift-core state. The magnetic moment is ordered antiferromagnetically along $`z`$ axis, and the amplitude is almost uniform taking its maximum at each site where the moment vanishes to form the antiferromagnetic state realized in the uniform state.

Stable structures near $`\mathrm{\Omega }_{c2}`$.— Having seen that the conventional Abrikosov lattice is never favored in the whole antiferromagnetic domain $`g_\mathrm{s}0`$, we now enumerate all the stable structures found near $`\mathrm{\Omega }_{c2}`$ to clarify their essential features. This rapidly rotating domain is especially interesting, because the same free energy also becomes relevant near $`\mathrm{\Omega }_{c2}`$ at $`T=0`$, as shown by Ho using a “mean-field quantum Hall regime” Ho01 . Thus, the conclusions obtained here are also applicable to the region at $`T=0`$.

Figure 2 displays the lowest free energy per unit volume as a function of $`g_\mathrm{s}`$ for $`\mathrm{\Omega }=0.95\mathrm{\Omega }_{c2}`$. The value of each $`n_\kappa `$ denotes the number of circulation quanta per unit cell. A special feature to be noted is that these various structures are energetically quite close to each other; for example, the $`n_\kappa =8`$ state at $`g_\mathrm{s}=0.02`$ is favored over the $`n_\kappa =3`$ state by a relative free-energy difference of order $`10^6`$. This fact suggests that we may realize quite a variety of metastable structures by a small change of the boundary conditions, the rotation process, etc.

Details of these structures are summarized as follows.

The $`n_\kappa =2`$ state of Eqs. (8) and (9) is stable for $`g_\mathrm{s}0.0671`$. We have already seen its basic features above.

The $`n_\kappa =3,5,7`$ states can be expressed compactly as

$`\mathrm{\Psi }_0={\displaystyle \underset{N}{}}{\displaystyle \underset{\nu =1}{\overset{\frac{n_\kappa 1}{2}}{}}}c_{N𝐪_\nu }^{(0)}\left(\psi _{N𝐪_\nu }\mathrm{e}^{2i\frac{\nu }{n_\kappa }\pi }\psi _{N𝐪_{n_\kappa \nu }}\right),`$ (10)
$`\mathrm{\Psi }_{\pm 1}={\displaystyle \underset{N}{}}[{\displaystyle \underset{\nu =1}{\overset{\frac{n_\kappa 1}{2}}{}}}c_{N𝐪_\nu }^{(\pm 1)}(\psi _{N𝐪_\nu }+\mathrm{e}^{2i\frac{\nu }{n_\kappa }\pi }\psi _{N𝐪_{n_\kappa \nu }})`$
$`+c_{N𝐪_{n_\kappa }}^{(\pm 1)}\psi _{N𝐪_{n_\kappa }}],`$ (11)
where $`N`$’s are even and $`𝐪_\nu =\frac{\nu }{n_\kappa }𝐛_1`$. These $`n_\kappa =3,5,7`$ states are stable for $`0g_\mathrm{s}0.0196`$, $`0.0313g_\mathrm{s}0.0613`$, and $`0.0613g_\mathrm{s}0.0671`$, respectively. Unlike the two component system considered by Mueller and Ho Ho02  where each component is specified by a single-$`𝐪`$ basis function, $`\mathrm{\Psi }_m`$ here is made up of multiple basis functions $`\psi _{N𝐪_\nu }`$ whose cores are shifted from each other by $`(\nu /n_\kappa )𝐚_2`$. Figure 3 displays the basic features of the $`n_\kappa =3`$ state. The lattice is hexagonal at $`g_\mathrm{s}=0`$, but deforms rapidly as $`g_\mathrm{s}`$ increases. The order-parameter amplitude is again almost constant, and the magnetic moment $`𝐌`$ has a complicated structure. These features are common to all the $`n_\kappa 3`$ states discussed here, although no details are presented for the other states. The lattice parameters $`(\beta ,\rho )`$ for $`n_\kappa =5,7`$ states are $`(0.1311\pi ,2.056)`$ and $`(0.1830\pi ,2.609)`$ for $`g_\mathrm{s}=0.05`$ and $`0.066`$, respectively, which change little in each relevant range of stability.

The remaining $`n_\kappa =8`$ state, stable over $`0.0196g_\mathrm{s}0.0313`$, is given by

$`\mathrm{\Psi }_0={\displaystyle \underset{N}{}}c_{N𝐪_1}^{(0)}\left(\psi _{N𝐪_1}+\psi _{N𝐪_2}\psi _{N𝐪_3}+i\psi _{N𝐪_4}\right),`$ (12)
$`\mathrm{\Psi }_{\pm 1}={\displaystyle \underset{N}{}}{\displaystyle \underset{\nu =5}{\overset{8}{}}}c_{N𝐪_\nu }^{(\pm 1)}\psi _{N𝐪_\nu },`$ (13)
where $`𝐪_1=\frac{𝐛_1+𝐛_2}{2}`$, $`𝐪_2=\mathrm{𝟎}`$, $`𝐪_3=\frac{𝐛_1𝐛_2}{4}`$, $`𝐪_4=\frac{𝐛_1+𝐛_2}{4}`$, $`𝐪_5=\frac{𝐛_1}{2}`$, $`𝐪_6=\frac{𝐛_2}{2}`$, $`𝐪_7=\frac{𝐛_1+𝐛_2}{4}`$, and $`𝐪_8=\frac{𝐛_1𝐛_2}{4}`$. The parameters $`(\beta ,\rho )`$ at $`g_\mathrm{s}=0.025`$ are $`(0.3317\pi ,1.049)`$, and changes only slightly in the above range of $`g_\mathrm{s}`$.

Concluding remarks.—We have performed extensive calculations on antiferromagnetic $`F=1`$ spinor vortices in rapid rotation. The conventional Abrikosov lattice is shown never favored. Each stable structure has almost constant order-parameter amplitude and a complicated magnetic-moment configuration, as shown in Figs. 1 and 3, for example. This means that any optical experiments to detect vortices by amplitude reduction will not be suitable for the spinor vortices. Instead, techniques to directly capture spatial magnetic-moment configurations will be required. The system has many metastable structures which are different in the number of circulation quanta per unit cell $`n_\kappa `$, but are quite close to each other energetically. Thus, domains to separate different structures may be produced easily. This degenerate feature is also present within the $`(\beta ,\rho )`$ space of a fixed $`n_\kappa `$, where $`\beta `$ is the vortex-lattice apex angle and $`\rho `$ is the length ratio of the two basic vectors. Put it differently, we can deform a stable lattice structure with a tiny cost of energy. These facts indicate that the spinor BEC’s can be a rich source of novel vortices realized with a small change in external parameters.

This research is supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology of Japan.