text:
# Renormalized Jellium model for charge-stabilized colloidal suspensions

## Abstract

We introduce a renormalized Jellium model to calculate the equation of state for charged colloidal suspensions. An almost perfect agreement with Monte Carlo simulations is found. Our self-consistent approach naturally allows to define the effective charge of particles at finite colloidal density. Although this quantity may differ significantly from its counterpart obtained from the standard Poisson-Boltzmann cell approach, the osmotic pressures for both models are in good agreement. We argue that by construction, the effective charge obtained using the Jellium approximation is more appropriate to the study of colloidal interactions. We also discuss a possibility of a fluid-fluid critical point and show how the new equation of state can be used to shed light on the surprising results found in recent sedimentation experiments.

In spite of the great effort invested in trying to understand the phase stability of colloidal suspensions, our knowledge of these complex systems is still quite rudimentary. It is curious to compare this situation with an earlier debate concerning the nature, or even the possibility of the liquid-gas phase separation in symmetric electrolytes. Now this debate is almost over, and the phase structure of a symmetric 1:1 electrolyte is well elucidated, although the universality class of the critical point is still being discussed. The Coulombic criticality —to distinguish it from the solvophobic criticality— cannot be observed in water and organic solvents of low dielectric permittivity must be used Fis94 . The phase separation results from the strong electrostatic correlations between the cations and anions of electrolyte Lev02 .

The phase stability of charged colloidal suspensions, on the other hand, is far from being well understood vRoij ; Linse ; Warren ; Diehl ; Deserno ; Trizac ; Tamashiro ; Dufreche . The strong charge and size asymmetry between the macroions and the microions present in a suspension makes it very difficult to apply to these systems the traditional tools of the liquid state theory. Nevertheless, in view of results on symmetric 1:1 electrolytes, it would not be very surprising if colloidal suspensions inside a solvent of sufficiently low dielectric permittivity $`ϵ`$ also presented a gas-liquid phase transition. Indeed such an instability has been observed in recent Monte Carlo simulations, see for example Lobaskin ; Linse  and references therein. What is much more surprising is that there are some experimental indications of an instability even in aqueous suspensions containing only monovalent counterions. Theoretical estimates of the strength of electrostatic correlations, for aqueous suspensions suggest that they should be too “hot” for an instability to set in. Nevertheless the experimental situation remains unclear Bla .

A number of theories have been proposed to address this unsettling experimental situation. A major drawback of most of these approaches is that they rely on uncontrolled approximations which have not been fully tested. However, numerical solution of the full non-linear Poisson-Boltzmann (PB) equation inside a spherical Wigner-Seitz (WS) cell Marcus  finds no indication of any thermodynamic instability Deserno . Of course, one can rightly question the reliability of the Wigner-Seitz cell PB model for the study of a fluid phase of a highly disordered suspension. On the one hand, since the main contribution to the osmotic pressure inside an aqueous suspension comes from the polyion-microions interactions, the cell might not be such a bad approximation. On the other hand, the cell fails to properly account for the colloid-colloid correlation, but these might not be of much importance for aqueous suspensions with monovalent counterions. Nevertheless, while the cell model is a good approximation for dense colloids, one should be very careful in extrapolating its findings to highly dilute suspensions. Clearly there is an urgent need for an accurate theory which will not rely on the cell approximation and which would be relevant for the study of colloidal phase stability. In this Letter we shall present such a theory. Our approach is similar to the Jellium approximation much used in the solid state physics.

Consider an aqueous suspension of colloidal particles of charge $`Zq`$ and radius $`a`$ in contact with a reservoir of monovalent salt at concentration $`c_s`$ and electrostatic potential $`\varphi _r=0`$ ($`q`$ is the elementary charge). The number of counterions and coions inside the suspension is determined by the thermodynamic equilibrium. While the colloidal particles are more or less uniformly distributed throughout the solution —we are mainly interested in the small density regime—, the positions of counterions and coions are strongly correlated with the positions of polyions. As a leading order approximation we can, therefore, take the polyion-polyion correlation function to be $`g_{pp}=1`$ Beresford  while the exact polyion-counterion and polyion-coion correlation functions are $`g_{p\pm }=e^{\beta w_\pm (r)}`$, where $`w_\pm (r)`$ are the polyion-counterion and the polyion-coion potentials of mean force and $`\beta =1/(k_BT)`$ is the inverse temperature.

Choosing the coordinate system in such a way that it is centered on top of one of the colloidal particles, the electrostatic potential satisfies the Poisson equation

$$^2\varphi =\frac{4\pi }{ϵ}\rho _q(r).$$
(1)
The charge density is $`\rho _q(r)=\rho _{\text{back}}+q\rho _+(r)q\rho _{}(r)`$, where, $`\rho _{\text{back}}=Z_{\text{back}}q\rho _p`$, and $`\rho _p`$ is the mean density of colloids inside suspension. The background charge is excluded from the colloidal interior. Naively one can suppose that $`Z_{\text{back}}=Z`$. This, however, is not correct and the bare charge must be renormalized in such a way as to lead to a self-consistent solution of Eq. (1), as discussed below.

Approximating the potential of mean force by the electrostatic potential, the local concentration of counterions and coions inside a suspension is,

$$\rho _\pm (r)=c_se^{\beta q\varphi (r)}.$$
(2)
It is important to keep in mind that in order for suspension to be neutral, the electrostatic potential at infinity (bulk) cannot vanish but must saturate to a value $`\varphi (\mathrm{})=\varphi _D`$ given by $`c_s\mathrm{sinh}[\beta q\varphi _D]=\rho _pZ_{\text{back}}`$. There exists, therefore, an electrostatic potential difference between the suspension and the salt reservoir. In the biophysics literature, this potential difference is known as the Donnan potential, and is partially responsible for the biological cell trans-membrane potential.

Far away from the colloidal surface, the electrostatic potential reduces to the familiar Debye-Hückel expression

$$\varphi (r)=\varphi _D\frac{Z_{\text{eff}}q}{ϵ(1+\kappa a)r}e^{\kappa (ra)},$$
(3)
where $`\kappa ^2=4\pi \lambda _B[\rho _+(\mathrm{})+\rho _{}(\mathrm{})]`$ and $`\lambda _B=\beta q^2/ϵ`$ denotes the Bjerrum length. The value of $`Z_{\text{eff}}`$ is determined self-consistently from the numerical solution of Eq. (1) so that $`Z_{\text{eff}}(Z,Z_{\text{back}},c_s)=Z_{\text{back}}`$. This renormalization of background charge is a consequence of counterion condensation. It is important to keep in mind that for highly charged colloidal particles, $`Z_{\text{eff}}`$ is not equal to the bare colloidal charge. Furthermore, using the contact theorem and the vanishing of the electric field as $`r\mathrm{}`$, the osmotic pressure within the suspension takes a simple form

$$\beta P=\rho _p+\rho _+(\mathrm{})+\rho _{}(\mathrm{})2c_s=\rho _p+\sqrt{Z_{\text{eff}}^2\rho _p^2+4c_s^2}2c_s.$$
(4)
It is a nice feature of the Jellium approximation that once the effective charge is determined, the Debye length $`1/\kappa `$ and the osmotic pressure both follow directly.

In Fig. 1 we compare the osmotic pressures calculated using the renormalized Jellium approximation to the ones obtained within the WS cell (for both approaches, the effective charges and pressures become independent of $`Z`$ when the latter is large enough, approximately $`Z>20a/\lambda _B`$; Fig. 1 has been obtained under this condition of saturation, which is usually met in colloidal suspensions). A surprisingly good agreement is found between the two theories, with a discrepancy only for volume fractions $`\eta >0.15`$. In Fig 2, the osmotic pressure calculated in the renormalized Jellium model is compared to the results of recent Monte Carlo simulations Linse  (where a model system of charged spherical macroions and point counterions interacting solely through hard sphere and Coulomb forces has been considered). The agreement is excellent, and justifies the neglect of colloid-colloid correlations in our treatment. Clearly, for high packing fractions (namely $`\eta >0.1`$) colloid-colloid correlations become important and should invalidate our approach.

In Fig. 3, the effective charges calculated using the cell model Alexander  and the renormalized Jellium theory are compared. At very low volume fractions ($`\eta <10^5`$) the saturated effective charge of a salt-free suspension is in perfect agreement with the value of saturated effective charge found within the WS cell and can be approximated by a simple equation,

$$Z_{\text{sat}}\frac{a}{\lambda _B}\left[\delta \gamma \mathrm{ln}(\eta )\right],$$
(5)
where $`\gamma 1`$ and $`\delta 2`$. For volume fractions $`\eta >10^4`$, there is a fairly strong disagreement between the effective charges predicted by the two theories, even if both models are in very good agreement for the value of the osmotic pressure.

If the effective charge is to be used to study the structural properties of colloidal suspensions (such as the structure factors or any other quantity requiring the knowledge of the effective interaction potential between two polyions), we argue that the Jellium effective charge is the more relevant quantity. The reason for this is that within the Jellium approximation, two colloids at large separations $`R`$ interact by the usual DLVO potential, as follows after some algebra from integrating the stress tensor over a colloid’s surface (see also Lev02 )

$$V(R)=\frac{Z_{\text{eff}}^2q^2}{ϵ(1+\kappa a)^2R}e^{\kappa (2aR)}.$$
(6)
This is not the case for two colloids inside a WS cell for which the interaction potential is a very complicated function of separation and can only be calculated numerically. In fact, for separations larger than the inter-particle distance, it is even difficult to properly define what one means by an interaction potential between colloids within the WS formalism, since the overall charge neutrality results in a vanishing electrostatic interaction between the two cells.

To assess the quantitative validity of our approach, we compare our (saturated) effective charges to those deduced in recent experiments using a torsional resonance spectroscopy Wette . For colloidal volume fractions in the interval $`10^3<\eta <3.10^2`$, $`Z_{\text{sat}}\lambda _B/a`$ was found to be close to $`6`$ Wette . This is in much better agreement with the renormalized Jellium model, which finds for these volume fractions $`Z_{\text{sat}}\lambda _B/a`$ varying between $`6.77.0`$, as compared to the WS cell prediction of $`8.79`$. We also note that a value close to 6 has been reported for this density regime in the Monte Carlo study presented in Ref Stevens .

We also emphasize that within the Jellium model, the inverse screening length $`\kappa `$ is naturally related to the effective salt density \[see the definition below Eq. (3)\]. This should be contrasted with the cell model, for which there is no simple connection between the two quantities.

The Jellium model predicts stability of a charged colloidal suspension against a fluid-fluid phase separation. This reinforces the cell picture, where no instability is found. Our approach, however, neglects microionic correlations, which become important at high electrostatic couplings Rouzina ; Schklovskii ; Lev02 , more precisely when $`\mathrm{\Gamma }=v^2\lambda _B/d`$ exceeds a threshold close to 2, $`d`$ being a characteristic distance between microions in the electric double-layer and $`v`$ their valency. In practice, however, ionic hydration puts a lower bound to $`d`$, which prevents the high coupling regime from ever being reached in water with monovalent counterions. Under these conditions our approach should, therefore, be quite reliable. Alternatively, when dealing with point particles, $`d`$ may be estimated as $`d(4\pi a^2v/Z)^{1/2}`$ Rouzina ; Lev02 . For the Monte Carlo simulations reported in Fig. 2, $`\mathrm{\Gamma }=2`$ therefore corresponds to an instability threshold $`\lambda _B/a1.1`$. Beyond this point, the microionic correlations destabilize the system and lead to a fluid-fluid separation Lobaskin ; Linse . Nevertheless, the pressures obtained within the renormalized Jellium approximation are in very good agreement with the Monte Carlo simulations, even at fairly large electrostatic couplings (in Fig. 2, we have included the highest coupling for which the pressure was computed in Linse ). We therefore expect that the Jellium approximation, suitably corrected to include the counterion induced attraction between the colloidal particles present at strong couplings, might be sufficient to account for the thermodynamic instability in suspensions containing multivalent ions. Work along these lines is in progress.

Finally, it is interesting to speculate how the renormalized Jellium model can help to shed new light on the problem of sedimenting colloidal dispersions Piazza ; JPH ; Simonin ; Lowen ; Tellez . Recently, Philipse and Koenderink (PK) PK  observed strongly inflated concentration profiles for charged mono-dispersed colloidal particles in absolute ethanol. Use of ethanol, instead of water, allows to produce highly deionized suspensions with “salt” concentrations on the order of $`10^9`$ M. The renormalized Jellium predicts that for this salt concentration, an infinitely dilute colloid would have an effective charge of $`Z_{\text{eff}}20a/\lambda _B`$. It is difficult to know what is precisely the bare charge of colloidal particles inside ethanol, however since the Bjerrum length in ethanol is $`\lambda _B2.3`$ nm, a fairly small bare charge on the order of a few hundred electrons is enough to place the particles (with radius $`a90`$nm PK ) in the saturation regime. Using this observation we can partially account for the observations of PK.

Static equilibrium of a suspension in a gravitational field requires that

$$\frac{dP}{dx}=mg\rho _p,$$
(7)
where $`m`$ is the colloidal mass (corrected for buoyancy), $`g`$ the gravitational acceleration, and $`x`$ the vertical displacement. It is convenient to define the gravitational length as $`l_g=1/(mg\beta )`$, which for the experiments of PK is $`l_g0.2`$mm PK . Substituting the expression for the osmotic pressure Eq.(4) into Eq.(7) and using the dependence of the saturated effective charge on the colloidal volume fraction, Eq.(7) can be integrated. We find that for colloidal volume fractions $`\eta >\eta _1=10^5`$ (i.e. in the counterion dominated regime):

$$\mathrm{ln}(\eta /\eta _0)\frac{1}{2}\left[\mathrm{ln}^2(\eta )\mathrm{ln}^2(\eta _0)\right]=\frac{\lambda _Bx}{al_g},$$
(8)
where $`\eta _0`$ is the reference volume fraction at the $`x`$ origin. This can be taken as $`\eta _00.01`$, the point at which the hard-core effects are completely negligible, and one can be sure to be looking at the dilute tail of the sedimenting profile. Note that for small volume fractions, $`\eta \mathrm{exp}(\sqrt{2\lambda _Bx/al_g})`$ unlike the simple exponential barometric law. For $`\eta <\eta _\times 5\times 10^8`$, the salt resulting from solvent dissociation dominates over counterions, and we recover the usual barometric law $`\eta (x)\mathrm{exp}(x/l_g).`$ Results concerning the crossover regime $`\eta _\times <\eta <\eta _1`$ will be published elsewhere. It is interesting to estimate the extent of the sedimentation profile predicted by the renormalized Jellium model. Since the concentrations of colloids on the order of $`\eta _m10^5`$ are still detectable, we find that the profile extends distance $`xal_g\mathrm{ln}^2(\eta _m)/(2\lambda _B)2600l_g`$. For the colloids used by PK, this is almost $`50`$ cm! Hence, the observed inflation of the sedimentation profile.

To conclude, we have proposed a self-consistent renormalized non-linear Jellium model that constitutes an alternative to the widely used PB cell approach for charged colloidal suspensions. Surprisingly, the equations of state within the two theories turn out to be very close, even though that the approximations involved to account for finite colloidal density are very different (finite cell against a renormalized background). Our results point to the relevance of the cell picture even at extremely low densities. This reinforces the argument that a fluid-fluid phase instability is impossible, for highly charged colloids in water at room temperature, as long as the microions are monovalent \[which allows to neglect microionic correlations, thereby identifying the potential of mean force with the electrostatic potential, see Eq. (2)\]. Our approach also allows to define in a natural way, not only the effective charge of the macroions, but also the effective screening length; these quantities are experimentally and conceptually more relevant than those obtained within the cell approach. Finally, our equation of state for deionized systems helps in understanding recent sedimentation experiments where “anomalous” density profiles have been reported.

This work was supported in part by the Brazilian agencies CNPq and FAPERGS and by the french CNRS.


format: text

text:
# The Formation and Early Evolution of Protostellar Accretion Disks

## 1. Introduction

Protostellar disks are of much interest in the study of star formation since it is likely that most of the mass assembled in a typical low-mass young stellar object (YSO) is accreted through a disk. Furthermore, such disks are the incubators of planetary systems, so their properties are also directly relevant to the process of planet formation. Circumstellar disks have been detected in $`2550\%`$ of pre–main-sequence stars in nearby dark clouds (e.g., Beckwith & Sargent 1993). At the time when the protostars become visible as classical T Tauri stars, their disk masses (as inferred from observations of dust emission and from spectroscopic measurements) are usually $`<10\%`$ of the central mass.

Rotationally supported circumstellar disks evidently originate in the collapse of self-gravitating, rotating, molecular cloud cores. Molecular line observations (e.g., Goodman et al. 1993; Kane & Clemens 1997) have established that a majority of dense ($`>10^4\mathrm{cm}^3`$) cloud cores show evidence of rotation, with angular velocities $`3\times 10^{15}10^{13}\mathrm{s}^1`$ that tend to be uniform on scales of $`0.1\mathrm{pc}`$, and with specific angular momenta in the range $`4\times 10^{20}3\times 10^{22}\mathrm{cm}^2\mathrm{s}^1`$. The cores can transfer angular momentum to the ambient gas through torsional Alfvén waves (by the process of magnetic braking), and this mechanism also acts to align their angular momentum vectors with the local large-scale magnetic field (e.g., Mouschovias & Ciolek 1999). This alignment can occur on the dynamical timescale and hence could be achieved even in cores whose lifetimes are of that order (as in the hierarchical-ISM scenario of Elmegreen 2000). Once dynamical collapse is initiated and a core goes into a near–free-fall state, the specific angular momentum is expected to be approximately conserved, resulting in a progressive increase in the centrifugal force that eventually halts the collapse and gives rise to a rotationally supported disk on scales $`10^2\mathrm{AU}`$. These expectations are consistent with the results of molecular-line interferometric observations of contracting cloud cores (e.g., Ohashi et al. 1997; Belloche et al. 2002).

This contribution summarizes recent advances in semianalytic modeling of the formation and early evolution of protostellar accretion disks. It has now become possible to study these processes in the context of a realistic scenario of star formation in molecular cloud cores that are threaded by a dynamically significant magnetic field. A more detailed account of this work is given in Krasnopolsky & Königl (2002; hereafter KK02).

## 2. Modeling Framework

Numerical simulations of magnetically supported clouds have demonstrated that the gas rapidly contracts along the field lines and maintains force equilibrium along the field even during the collapse phase (e.g., Fiedler & Mouschovias 1993; Galli & Shu 1993), including in cases where the clouds are initially elongated in the field direction (e.g., Nakamura, Hanawa, & Nakano 1995; Tomisaka 1996). This motivates treating the collapse as being quasi one-dimensional.111The complementary approach of axisymmetric, 2D collapse simulations in which new mass is added to the system only from above and below the disk plane tends to produce disk-to-YSO mass ratios $`1`$, much higher than typically observed.

To obtain semianalytic solutions, KK02 adopted the assumption of self-similarity in space and time, with a similarity variable $`xr/Ct`$ (where $`r`$ is the distance from the origin, $`C`$ is the isothermal speed of sound, and $`t`$ is the time).222For a typical sound speed $`C=0.19\mathrm{km}\mathrm{s}^1`$, $`x=1\{400,4000\}\mathrm{AU}`$ at $`t=\{10^4,10^5\}\mathrm{yr}`$. This assumption is motivated by the fact that core collapse is a multiscale problem, which is expected to assume a self-similar form away from the outer and inner boundaries and not too close to the onset time (e.g., Shu 1977; Hunter 1977). This behavior has been verified by previous numerical and semianalytic treatments of restricted core-collapse problems – with/without rotation and with/without magnetic fields. The assumption of isothermality, which underlies the $`C=const`$ ansatz, is justified mainly by the fact that thermal stresses do not play a major role in the dynamics of the collapsing core.

Molecular cloud cores are known to be weakly ionized. Therefore, even though they also carried out reference calculations under the assumption of ideal MHD, KK02 incorporated ambipolar diffusion into the model: the magnetic field lines are frozen into the charged particle component (ions, electrons, grains) and couple to the dominant neutral component through ion–neutral drift. Although the drift velocity is negligible during the early phase of the dynamical collapse, ambipolar diffusion becomes important within the gravitational “sphere of influence” of the YSO once the central mass begins to grow (Ciolek & Königl 1998; Contopoulos, Ciolek, & Königl 1998). When the incoming matter enters this region, it decouples from the field and continues moving inward. The decoupling front, in turn, moves outward and steepens into a C-type ambipolar diffusion (AD) shock (the existence of which was first predicted by Li & McKee 1996).

To incorporate ambipolar diffusion into the self-similarity formulation, it is necessary to assume that the ion density scales as the square root of the neutral density: $`\rho _i=K\rho ^{1/2}`$. As discussed in KK02, this is a good approximation for the core-collapse problem: it applies on both ends of a density range spanning $`8`$ orders of magnitude, which applies roughly on radial scales $`1010^4\mathrm{AU}`$, with $`K`$ varying by only 1 order of magnitude across this interval.

The transition from a nearly freely falling, collapsing core to a quasi-stationary, rotationally supported disk involves a strong deceleration in acentrifugal shock. This shock is distinct from the ambipolar-diffusion shock mentioned above: it typically occurs at a different radius and is hydrodynamic, rather than hydromagnetic, in nature.

To allow mass to accumulate at the center in a 1D, rotating-core collapse, an angular momentum transport mechanism must be present. In their basic model, KK02 assumed that vertical transport through magnetic braking continues to operate also during the collapse phase of the core evolution. To incorporate this mechanism into the self-similar model, it is necessary to assume that $`V_{\mathrm{A},\mathrm{ext}}`$, the Alfvén speed in the external medium, is a constant.333A nearly constant value $`V_{\mathrm{A},\mathrm{ext}}1\mathrm{km}\mathrm{s}^1`$ is, in fact, indicated in molecular clouds in the density range $`10^310^7\mathrm{cm}^3`$ (e.g., Crutcher 1999). KK02 verified that, in their derived solutions, magnetic braking indeed dominates the most likely alternative angular-momentum transport mechanisms — magnetorotational instability-induced turbulence and gravitational torques. However, they also found that angular momentum transport by a centrifugally driven magnetic disk wind arises naturally (and may dominate) in their fiducial disk solutions. They went on to show that the latter mechanism may be incorporated into the model without significantly modifying the basic formulation.

## 3. Self-Similar Model

KK02 employed a set of vertically-integrated thin-disk equations. In their formulation (using cylindrical coordinates $`r,\varphi ,z`$), they took the vertical magnetic field component to be constant (except when evaluating $`B_z/z`$) and assumed that $`B_r`$ and $`B_\varphi `$ increase $`z`$ in the disk; all terms $`𝒪(H/r)`$ (where $`H`$ is the disk scale height) were neglected except in the azimuthal current-density term $`[B_{r,s}H(B_z/r)]`$ (where the subscript $`s`$ denotes a surface value). Furthermore, a monopole approximation was employed for the radial gravity $`g_r`$ and for $`B_{r,s}`$ (e.g., Li & Shu 1997): $`g_r=GM(r,t)/r^2`$, $`B_{r,s}=\mathrm{\Psi }(r,t)/2\pi r^2`$, relating these quantities to, respectively, the mass and magnetic flux enclosed within $`r`$.

In the self-similar formulation, the various physical quantities are expressed as dimensionless functions of the similarity variable $`x`$ in the following fashion:

$$H(r,t)=Cth(x),\mathrm{\Sigma }(r,t)=(C/2\pi Gt)\sigma (x),$$

$$V_r(r,t)=Cu(x),V_\varphi (r,t)=Cv(x),$$

$$g_r(r,t)=(C/t)g(x),J(r,t)=C^2tj(x),$$

$$M(r,t)=(C^3t/G)m(x),\dot{M}(r,t)=(C^3/G)\dot{m}(x),$$

$$𝐁(r,t)=(C/G^{1/2}t)𝐛(x),\mathrm{\Psi }(r,t)=(2\pi C^3t/G^{1/2})\psi (x),$$

where $`\mathrm{\Sigma }`$ is the surface mass density, $`𝐕`$ is the velocity, $`J`$ is the specific angular momentum, and $`\dot{M}`$ is the mass accretion rate. From the assumption of vertical hydrostatic equilibrium, one can derive a quadratic equation for $`h(x)`$, whose solution is

$$h=\frac{\widehat{\sigma }x^3}{2\widehat{m}_c}\left[1+\left(1+\frac{8\widehat{m}_c}{x^3\widehat{\sigma }^2}\right)^{1/2}\right],$$

where $`\widehat{m}_cm_cx^3b_{r,s}(db_z/dx)/\sigma `$ (with $`m_c`$ being the central-mass eigenvalue) and $`\widehat{\sigma }\sigma +(b_{r,s}^2+b_{\varphi ,s}^2)/\sigma `$.

The adopted initial conditions were $`\sigma \frac{A}{x},b_z\sigma /\mu _0,uu_0,vv_0`$ as $`x\mathrm{}`$, with the parameter values $`A=3`$, $`\mu _0=2.9`$, $`u_0=1`$ based on numerical simulations of nonrotating magnetic cores and with the choice of the value of $`v_0`$ motivated by the measured angular velocities in molecular cloud cores.

The behavior in the limit $`x0`$ (corresponding to $`r0`$ at a fixed $`t`$) can be derived from the constituent equations. In particular, the asymptotic behavior of an ambipolar diffusion-dominated circumstellar disk is given by

$`\dot{m}`$ $`=`$ $`m=m_c,`$
$`j`$ $`=`$ $`m_c^{1/2}x^{1/2},`$
$`u`$ $`=`$ $`w=(m_c/\sigma _1)x^{1/2},`$
$`\sigma `$ $`=`$ $`{\displaystyle \frac{(2\eta /3\delta )(2m_c)^{1/2}}{[1+(2\eta /3\delta )^2]^{1/2}}}x^{3/2}`$
$``$ $`\sigma _1x^{3/2},`$
$`b_z`$ $`=`$ $`b_{\varphi ,s}/\delta =[m_{c}^{}{}_{}{}^{3/4}/(2\delta )^{1/2}]x^{5/4},`$
$`b_{r,s}`$ $`=`$ $`\psi /x^2=(4/3)b_z,`$
$`h`$ $`=`$ $`\{2/[1+(2\eta /3\delta )^2]m_c\}^{1/2}x^{3/2}.`$
There are 4 model parameters that can be varied to explore the solution space: $`\delta |B_{\varphi ,s}|/B_z`$ (the adopted cap on the azimuthal field strength), $`\alpha C/V_{\mathrm{A},\mathrm{ext}}`$, $`v_0V_{\varphi ,0}/C`$, and $`\eta \tau _{\mathrm{ni}}(4\pi G\rho )^{1/2}`$ (where $`\tau _{ni}1/\rho _i`$ is the neutral–ion momentum-exchange time).

## 4. Results

### 4.1. Fiducial Solution

This solution corresponds to $`\eta =1`$, $`v_0=0.73`$, $`\alpha =0.08`$, and $`\delta =1`$ (which yield $`m_c=4.7`$). In this case the initial rotation is not very fast and the braking is moderate, leading to the formation of a disk (with outer boundary at the centrifugal shock radius $`x_c=1.3\times 10^2`$) within the ambipolar-diffusion region (enclosed by the AD shock radius $`x_a=0.4130x_c`$). One can distinguish the following main flow regimes (see Fig. 1):

* Outer region ($`x>x_a`$): ideal-MHD infall.
* AD shock: resolved as a continuous transition (but may in some cases contain a viscous subshock); KK02 estimated $`x_a\sqrt{2}\eta /\mu _0.`$
* Ambipolar diffusion-dominated infall ($`x_c<x<x_a`$): near free-fall controlled by the central YSO’s gravity.
* Centrifugal shock: its location depends sensitively on the diffusivity parameter $`\eta `$, which affects the amount of magnetic braking for $`x<x_a`$; KK02 estimated $`x_c(m_cv_0^2/A^2)\mathrm{exp}[(2^{3/2}m_c/\mu _0)^{1/2}\eta ^{3/2}]`$.
* Keplerian disk ($`x<x_c`$): asymptotic behavior is approached after a transition zone representing a massive ring (of width $`0.1x_c`$ and mass $`8\%`$ of the disk mass within $`x_c`$, which in turn is $`<5\%`$ of $`m_c`$).

The asymptotic $`x0`$ solution (see § 3) implies that the angle between the meridional projection of $`𝐁`$ and the rotation axis is radially constant and equal to $`53^{}`$, which exceeds the minimum value of $`30^{}`$ for launching a centrifugally driven wind from a “cold” Keplerian disk (Blandford & Payne 1982). This feature of the solution is attractive in view of the fact that centrifugally driven disk winds are a leading candidate for the origin of the bipolar outflows that are frequently observed to emanate from YSOs (e.g., Königl & Pudritz 2000).

### 4.2. Limiting Cases: Fast Rotation and Strong Braking

By modifying the model parameters, one can study the range of possible behaviors in collapsing cores. Figure 2 shows two limiting cases, which bracket the fiducial solution.

The fast rotation case differs from the fiducial solution primarily in having a large initial-rotation parameter ($`v_0=1.5`$). It has the following distinguishing features:

* The centrifugal shock is located within the self-gravity–dominated (and ideal-MHD) region; a back-flowing region is present just behind the shock.
* The central mass is comparatively small ($`m_c=0.5`$), giving rise to a non-Keplerian outer disk region.
* The ideal-MHD/ambipolar-diffusion transition occurs behind the centrifugal shock and is gradual rather than sharp.

The strong braking case ($`\eta =0.5`$, $`v_0=1`$, $`\alpha =10`$, and $`\delta =10`$; yielding $`m_c=5.9`$) is characterized by large values of the braking parameters $`\alpha `$ and $`\delta `$. It is distinguished by having

* no centrifugal shock (or circumstellar disk); the $`x0`$ behavior resembles that of the nonrotating collapse solution of Contopoulos et al. (1998).

## 5. Some Implications of the Model

The formation of rotationally supported circumstellar accretion disks basically resolves the angular momentum problem in star formation (although the exact value of the YSO angular momentum is determined by processes near the stellar surface that are not included in this model). In particular, the derived solutions demonstrate that angular momentum transport can be sufficiently efficient to allow most of the inflowing mass to end up (with effectively no angular momentum) at the center, with the central mass dominating the dynamics well beyond the outer edge of the disk even as the inflow is still in progress. These solutions reveal that the ambipolar-diffusion shock, even though it is usually located well outside the region where the centrifugal force becomes important, helps to enhance the efficiency of angular momentum transport through the magnetic field amplification that it induces. The revitalization of ambipolar diffusion behind the AD shock in turn goes a long way toward resolving the magnetic flux problem in star formation (as already pointed out by Ciolek & Königl 1998 and Contopoulos et al. 1998).

To the extent that self-similarity is a good approximation to the situation in collapsing cloud cores, it is conceivable that T Tauri (Class II) protostellar systems, whose disk masses are typically inferred to be $`<10\%`$ of the central mass, have had a similarly low disk-to-star mass ratio also during their earlier (Class-0 and Class-I) evolutionary phases. It would be interesting to test this possibility by observations. The model also predicts that, in cases where magnetic braking is particularly strong, essentially all the angular momentum is removed well before the inflowing gas reaches the center. Such systems may correspond to slowly rotating YSOs that show no evidence of a circumstellar disk (e.g., Stassun et al. 1999; 2001). Another distinguishing characteristic of the solutions is the appearance of an ambipolar-diffusion and/or a centrifugal shock. The implied processing of the disk material in these shocks (particularly the latter one) may have implications to the composition of protoplanetary disks (e.g., the annealing of silicate dust; see Harker & Desch 2002).

The diffusive Keplerian disk models are by and large magnetorotationally stable, basically because the matter/field coupling is generally too weak to allow the instability to grow. (The well-coupled surface layers should typically also be stable because of strong magnetic squeezing; e.g., Wardle & Königl 1993.) Furthermore, the Toomre stability criterion to fragmentation ($`Q_{\mathrm{Toomre}}>1`$) is well satisfied for the rotationally supported disk solutions (except in the outer layers of fast-rotation models). Angular momentum transport by gravitational torques is unlikely to be important under these conditions (e.g., Lin & Pringle 1987). However, as noted in § 2 (see also § 4.1), angular momentum transport by a centrifugally driven wind may play a key role. KK02 found that the steady-state, radially self-similar disk-wind solution of Blandford & Payne (1982) can be naturally incorporated into the asymptotic ambipolar-diffusion disk solution given in § 3, making it possible to study the effects of wind angular-momentum and mass removal from the disk and to better constrain the relevant parameters of a combined disk/wind model. In a preliminary analysis, they inferred that the asymptotic solution evidently corresponds to the weakly coupled disk/wind configurations discussed by Li (1996).

#### Acknowledgments.

My work on star formation was initiated in 1980 when I came to Berkeley as a postdoc and started interacting with Chris, Dave, and Frank. Much of what I have learned about star formation in the ISM during the subsequent years has derived from our honorees’ work and was inspired by their physical insights and unique styles. I wish them all many years of continued productive leadership in this field. The research reported in this contribution was supported in part by NASA grant NAG5-3687.

## References

Beckwith, S. V. W., & Sargent, A. I. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine (Tucson: Univ. Arizona Press), 521

Belloche, A., André, P., Despois, D., & Blinder, S. 2002, A&A, 393, 927

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Ciolek, G. E., & Königl, A. 1998, ApJ, 504, 257

Contopoulos, I., Ciolek, G. E., & Königl, A. 1998, ApJ, 504, 247

Crutcher, R. M. 1999, ApJ, 520, 706

Elmegreen, B. G. 2000, ApJ, 530, 277

Fiedler, R. A., & Mouschovias, T. Ch. 1993, ApJ, 415, 680

Galli, D., & Shu, F. H. 1993, ApJ, 417, 243

Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528

Harker, D. E., & Desch, S. J. 2002, ApJ, 565, L109

Hunter, C. 1977, ApJ, 218, 834

Kane, B. D., & Clemens, D. P. 1997, AJ, 113, 1799

Königl, A., & Pudritz, R. E. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell (Tucson: Univ. Arizona Press), 759

Krasnopolsky, R., & Königl, A. 2002, ApJ, 580, 987 (KK02)

Li, Z.-Y. 1996, ApJ, 465, 855

Li, Z.-Y., & McKee, C. F. 1996, ApJ, 464, 373

Li, Z.-Y., & Shu, F. H. 1997, ApJ, 475, 237

Lin, D. N. C., & Pringle, J. E. 1987, MNRAS, 225, 607

Mouschovias, T. Ch., & Ciolek, G. E. 1999, in The Origins of Stars and Planetary Systems, ed. C. J. Lada & N. D. Kylafis (Dordrecht: Kluwer), 305

Nakamura, F., Hanawa, T., & Nakano, T. 1995, ApJ, 444, 770

Ohashi, N., Hayashi, M., Ho, P. T. P., Momose, M., Tamura, M., Hirano, N., & Sargent, A. I. 1997, ApJ, 488, 317

Shu, F. H. 1977, 214, 488

Stassun, K. G., Mathieu, R. D., Mazeh, T., & Vrba, F. J. 1999, AJ, 117, 2941

Stassun, K. G., Mathieu, R. D., Vrba, F. J., Mazeh, T., & Henden, A. 2001, AJ, 121, 1003

Tomisaka, K. 1996, PASJ, 48, 701

Wardle, M., & Königl, A. 1993, ApJ, 410, 218


text:
# An extragalactic HII region in the Virgo clusterBased on observations obtained with the Loiano telescope belonging to the University of Bologna (Italy) and with the Calar Alto observatory operated by the Centro Astronomico Hispano Aleman (Spain).

## 1 Introduction

The presence of a diffuse intra-cluster light component was for the first time proposed by Zwicky(zwicky (1951)), who detected a light excess between the galaxies in the Coma cluster. In the last twenty years, the presence of intra-cluster starlight (ICL) has been inferred from surface brightness measurements (Bernstein et al. bernstein (1995)) and from the detection of individual extragalactic stars (Ferguson et al. ferguson (1998)), planetary nebulae (Arnaboldi et al. arnaboldi96 (1996)), and HII regions (Gerhard et al. gerhard (2002)). The study of ICL and its origin could shed light on the evolution of galaxies in clusters, a complex process involving various mechanism such as tidal interactions (Gnedin gnedin (2003)), galaxy harassment (Moore et al. moore (1996)) and ram pressure (Gunn & Gott gunn (1972)). Gravitational interactions between cluster galaxies may strip gas and stars from these systems contributing to the ”pollution” of the intra-cluster medium (ICM) and forming, in some cases, extragalactic HII regions as the ones recently observed in the Abell cluster 1367 (Sakai et al. sakai (2002), Gavazzi et al. 2003a ). However the interstellar medium may also be removed from fast-moving galaxies by the ram pressure mechanism (Quilis et al. quilis (2000), Gavazzi et al. gavazzi01 (2001)). Some of this gas may collapse and form stars, which would contribute to the ICL. Moreover extragalactic massive stars might explode as supernovae enriching the ICM. Thus the ICL properties should provide important clues not only on the dynamical state of clusters and galaxies (Miller miller (1983)), but also on their chemical evolution. However, only one intra-cluster HII region candidate has been confirmed spectroscopically so far (Gerhard et al. gerhard (2002)).
In this work we present the result of spectroscopic observations of a second extragalactic HII region in the Virgo cluster, confirming its discovery by Cortese et al. (cortese (2003)).
We assume for the Virgo cluster a distance modulus of 31.15, corresponding to a distance of 17 Mpc.

## 2 Sample selection

With the aim of detecting faint emission line objects in the Virgo cluster region, we visually inspected the net ($`\mathrm{H}\alpha +[\mathrm{NII}]`$) frames obtained at the Isaac Newton Telescope by Boselli et al. (boselli (2002)). They observed nine different regions centered on BCD galaxies in the Virgo cluster, covering a total area of $``$2.5 square degrees with an $`\mathrm{H}\alpha `$ flux limit $`10^{15}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$. Since the seeing conditions during the observations varied field-to-field and even between ON and OFF-band observations, a conservative criterion was assumed for discriminating real from spurious detections, namely, we consider real those sources with $`EW>100\mathrm{\AA }`$ and with a $`S/N`$ ratio $`>`$ 5 in the net frame. Only six objects were found to meet this criterion. As shown in Table 3, they are all faint sources with detectable $`r^{}`$ continuum emission ($`20.1m_r^{}21.6`$) suggesting that they are not extragalactic planetary nebulae, since these types of sources have no continuum emission (Arnaboldi et al. arnaboldi03 (2003), Feldmeier et al. feldmeier (2003)). The emitting line systems detected in the narrow band filter could be extragalactic HII regions at the distance of the Virgo cluster, as the one recently discovered by Gerhard et al. (gerhard (2002)) or background emission-line galaxies. With the aim of understanding the real nature of these sources we undertook a spectroscopic follow-up.

## 3 Observations and data reduction

Long-slit, low dispersion spectra were obtained in several observing runs since 2002 using the imaging spectrograph BFOSC attached to the Cassini 1.5m telescope at Loiano (Italy), and with CAFOS attached to the 2.2m telescope of the Calar Alto Observatory (Spain). Table 1 lists the characteristics of the instrumentation in the adopted set-up. The journal of observations is presented in Table 2.
The observations at Loiano were performed using a 1.5 or 2.5 arcsec slit, depending on the seeing conditions, and the wavelength calibration was secured with exposures of HeAr lamps. The observations at Calar Alto were carried out with a 1.5 arcsec slit and the wavelength calibration was carried out with exposures of CdHe lamps. In all runs the observations were obtained in nearly photometric conditions or with thin cirrus.
The data reduction was performed in the IRAF111Image Reduction and Analysis Facility is written and supported by the National Optical Astronomy Observatories (NOAO) in Tucson, Arizona. NOAO is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under cooperative agreement with the National Science Foundation. environment. After bias subtraction, the cosmic rays were removed using the task COSMICRAYS and/or by visual inspection. The lamp wavelength calibration was checked against known sky lines. These were found within $`1\mathrm{\AA }`$ from their nominal positions, providing an estimate of the systematic uncertainty on the derived velocities of $`50\mathrm{km}\mathrm{s}^1`$. Sky subtraction was performed on the 2D spectra using the IRAF/LONGSLIT task BACKGROUND. BACKGROUND was used to fit the 2D spectra row-by-row with a first order Legendre function with $`3\sigma `$ rejection iteration. The sky was fitted in two windows of $`100\mathrm{pixels}`$, one on each side of the detected emission lines. After subtraction of sky background, one-dimensional spectra were extracted from the frames using $`APSUM`$. The apertures were limited to regions where the signal intensity was above 1 $`\sigma `$ of the sky noise. For the source *122603+130724* observations of the star FEIGE34 provided the flux calibration. The spectrophotometric standard was observed in nearly photometric conditions.

## 4 Results

The results of the measurements obtained in this work are listed in Table 3. The obtained spectra are shown in Fig 2. Five sources were found to be $`[\mathrm{OIII}]\lambda 5007`$ emitters at $`z0.31`$, only one source, *122603+130724*, belongs to the Virgo cluster. Cortese et al. (cortese (2003)) in their preliminary analysis, correctly interpreted the single emission line detected in the spectrum of *122603+130724* as due to $`\mathrm{H}\alpha `$ at the redshift of $`200\mathrm{km}\mathrm{s}^1`$, consistent with the velocity of VCC873 ($`232\mathrm{km}\mathrm{s}^1`$). Conversely the tentative interpretation that *122544+130740* might be another extragalactic HII region is ruled out by the present observations which show that the source is in fact a background galaxy and that the single emission line detected by Cortese et al. (cortese (2003)) is $`[\mathrm{OIII}]\lambda 5007`$ at $`z=0.315.`$
*122603+130724* lies at a projected distance of 3 kpc from VCC873 (see Fig.1). Boselli et al. (boselli (2002)) observed this bright galaxy in the narrow-band $`\mathrm{H}\alpha `$ filter and in the broad-band $`r^{}`$ and B filters. Spectrophotometric standards were observed only for $`\mathrm{H}\alpha `$ and B calibration. However from the narrow-band calibration we were able to calibrate the $`r^{}`$ frame too. The emitting line object *122603+130724* is clearly detected in both broad-band images as a point-like source with $`m_B22.25`$ and $`m_r^{}21.60`$. The uncertainties associated with these measures are $`0.25`$ mag due to contamination by VCC873. Not surprisingly the continuum is not detected in the spectrum of Fig. 2 because $`m_r^{}21.60`$ corresponds to a flux of $`f_r^{}\mathrm{6\; 10}^{18}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{\AA }^1`$, below the sensitivity limit $`f_r^{}10^{17}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2\mathrm{\AA }^1`$ of our observations.
From the narrow-band images we measured $`\mathrm{H}\alpha +[\mathrm{NII}]`$ flux $`=10^{14.56\pm 0.05}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$ and EW($`\mathrm{H}\alpha +[\mathrm{NII}]`$) $`700\pm 80\mathrm{\AA }`$. The uncertainty on the $`\mathrm{H}\alpha +[\mathrm{NII}]`$ flux includes three contributions: the Poisson photon counts error, the uncertainty on the background and the photometric uncertainty, which is assumed as 10% of the net flux. These errors were determined separately on the ON and OFF-band frames, and combined using the standard error propagation. The error in the equivalent width, is computed similarly to the flux uncertainty, except that the error on the absolute flux scale does not affect the equivalent width.
The flux-calibrated spectrum clearly shows three emission lines at the redshift of Virgo: $`\mathrm{H}\alpha `$, $`\mathrm{H}\beta `$ and $`[\mathrm{SII}]\lambda 6717,6731`$. The $`\mathrm{H}\alpha `$ line is blended with $`[\mathrm{NII}]\lambda 6583`$. Using the task SPLOT we deblended the two components providing a crude estimate of the line ratio $`[\mathrm{NII}]\lambda 6583/\mathrm{H}\alpha 0.17`$. The $`\mathrm{H}\alpha `$ flux obtained from the narrow-band images is consistent with the flux $`f_{H\alpha }10^{14.55}ergs^1cm^2`$ measured from the spectrum. Assuming a theoretical value of 2.86 for the Balmer decrement ($`T=10^4`$ K and $`n_e=100\mathrm{e}/\mathrm{cm}^3`$), we obtain $`A(H\beta )=1.67`$ and correct the line fluxes using the dereddening law of Cardelli et al. (cardelli (1989)). The dereddened $`\mathrm{H}\alpha `$ flux is $`f_{H\alpha }10^{14.10}ergs^1cm^2`$, translating into a luminosity $`L_{H\alpha }\mathrm{2.87\; 10}^{38}\mathrm{erg}\mathrm{s}^1`$. The $`[\mathrm{OII}]\lambda 3727,[\mathrm{OIII}]\lambda 5007\mathrm{and}4959`$ lines are undetected but we derive an upper limit to their intensity. The observed and reddening corrected line fluxes are listed in Table 4. We use the detected emission lines $`[\mathrm{NII}]\lambda 6583`$ and $`[\mathrm{SII}]\lambda 6717,6731`$ to obtain an estimate of the metallicity of *122603+130724* (Kewley & Dopita kewley (2002)). Unfortunately the ratio $`[\mathrm{NII}]\lambda 6583/[\mathrm{SII}]\lambda 6717,6731`$ is strongly dependent on the ionization parameter (see Fig. 4 in Kewley & Dopita kewley (2002)), thus we can only derive a lower limit $`8.612+\mathrm{log}(\mathrm{O}/\mathrm{H})`$, to the abundance of *122603+130724*.
The total number of H-ionizing photons estimated from (Osterbrock osterbrock (1989)) is:

$$Q(H^0)=\frac{\alpha _B}{\alpha _{H\alpha }^{eff}}\times \frac{L_{H\alpha }}{h\nu _{H\alpha }}2.8\times 10^{50}\mathrm{s}^1$$

where $`\alpha _\mathrm{B}/\alpha _{\mathrm{H}\alpha }^{\mathrm{eff}}2.96`$.
Thus the total mass and size of the gas cloud are:

$$M_{HII}=Q(H^0)\frac{m_p}{n_e\alpha _B}5150(\frac{n_e}{100\mathrm{cm}^3})^1\mathrm{M}_{}$$

and

$$r_{HII}=\left(\frac{M_{HII}(1+y^+)}{(4\pi /3)n_em_p}\right)^{1/3}8.2(\frac{n_e}{100\mathrm{cm}^3})^{2/3}\mathrm{pc}$$

where $`n_e`$ is the electron density, $`m_p`$ the proton mass and $`y^+0.1`$ is the fraction of ionized helium.
The total number of H-ionizing photons and the mass of the gas cloud are typical of giant HII regions (Stasinska stasinska (2002)). Moreover the detectable continuum emission of *122603+130724* and the low value of the ratio $`[\mathrm{OIII}]\lambda 5007/\mathrm{H}\alpha 0.6`$ rule out the possibility that this source is an intra-cluster planetary nebulae, confirming the hypothesis that it is an extragalactic HII region associated with the bright galaxy VCC873.
Recently Gerhard et al. (gerhard (2002)) found 17 candidate extragalactic HII regions in the Virgo cluster, among which only one, associated with VCC836, was confirmed spectroscopically. This source differs from the HII region presented in this work as it is smaller ($`r_{HII}3.5\mathrm{pc}`$), has a lower luminosity ($`L_{H\alpha }\mathrm{1.3\; 10}^{37}\mathrm{erg}\mathrm{s}^1`$) and metallicity ($`12+\mathrm{log}(\mathrm{O}/\mathrm{H})8.1`$) than *122603+130724*. Moreover the two giant galaxies to which the HII regions are associated have different properties. They are both edge-on star forming spiral galaxies but, while VCC836 is a Seyfert2 galaxy with a very extended emission-line region (Yoshida et al. yoshida (2002)), VCC873 is a non active object with no evidence of extended extra-planar $`\mathrm{H}\alpha `$ emission.
This evidence indicates the possibility that the two extragalactic HII regions could have different origins. Gerhard et al. (gerhard (2002)) suggest that the HII region associated with VCC836 might be a tidal debris of a past interaction with a gas-rich galaxy, as the very extended emission-line region associated with this bright galaxy indicates (Yoshida et al. yoshida (2002)). However Vollmer & Huchtmeier (vollmer (2003)) and Yoshida et al. (yoshida03 (2003)) propose an alternative interpretation in which the extended emission-line regions, and perhaps the extragalactic HII region, are formed by the interaction of VCC836 with the dense ICM of the Virgo cluster.
The mechanism responsible for the formation of the extragalactic HII region *122603+130724* is still not well understood, however the absence of a clear signature of interaction of VCC873 with another system rules out a tidal stripping scenario.
In recent years, deep HI observations of edge-on nearby spiral galaxies outlined the presence of a thick component (halo) of neutral hydrogen, with rotational velocity lower than the disk component, up to distances of 10-15 kpc from the disk (Swaters et al. swaters (1997)). Moreover large-scale radial inflow toward the galaxy center was revealed (Fraternali et al. fraternali (2001)). These observations suggest a complex gas circulation between the disk and the halo of star-forming spiral galaxies, supporting the galactic fountain model (e.g. Shapiro & Field shapiro (1976); Bregman bregman (1980)). This model predicts that the gas, ionized by supernova explosion and stellar winds, leaves the plane of the disk. After cooling, the ejected gas, presumably concentrated in clouds, falls back into the disk.
*122603+130724* lies within the halo of VCC873. Its velocity differs only by $`50\mathrm{km}\mathrm{s}^1`$ from the velocity observed in the disk of VCC873 near its location (Rubin et al. rubin (1999)), suggesting that *122603+130724* is gravitationally bound to the bright galaxy. The line ratio $`[\mathrm{NII}]\lambda 6583/[\mathrm{SII}]\lambda 6717,67311.2`$ of *122603+130724* is consistent with the metallicity obtained from the drift-scan mode spectrum of VCC873 itself, published by Gavazzi et al. (2003b , gavazzi04 (2004)), suggesting that this HII region formed from enriched material ejected from VCC873.
The proximity in position and velocity to the disk of VCC873 and the similar abundances observed in the two objects might indicate a galactic fountain origin for *122603+130724*. The expelled gas could have condensed into a cloud during its infall back to the disk, triggering star formation activity. If galactic fountains are as frequent as claimed by Fraternali et al. (fraternali (2001)), extragalactic HII regions possibly associated with this phenomenon might contribute significantly to the diffuse intra-cluster light.

###### Acknowledgements.

The TACS of the Loiano and Calar Alto telescopes are acknowledged for the generous time allocation to this project. This work could not have been completed without the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech under contract with NASA. We also made use of the GOLD Mine Database, operated by the Università degli Studi di Milano-Bicocca.


text:
# Nonlinear Band Gap Transmission in Optical Waveguide Arrays

## Abstract

The effect of nonlinear transmission in coupled optical waveguide arrays is theoretically investigated via numerical simulations on the corresponding model equations. The realistic experimental setup is suggested injecting the beam in a single boundary waveguide, linear refractive index of which ($`n_0`$) is larger than one ($`n`$) of other identical waveguides in the array. Particularly, the effect holds if $`\omega (n_0n)/c>2Q`$, where $`Q`$ is a linear coupling constant between array waveguides, $`\omega `$ is a carrier wave frequency and $`c`$ is a light velocity. Making numerical experiments in case of discrete nonlinear Schrödinger equation it is shown that the energy transfers from the boundary waveguide to the waveguide array above certain threshold intensity of the injected beam. This effect is explained by means of the creation and propagation of gap solitons in full analogy with the similar phenomenon of nonlinear supratransmission \[F. Geniet, J. Leon, PRL, 89, 134102, (2002)\] in case of discrete sine-Gordon lattice.

Nonlinear phenomena in large diversity of physical systems have a close relation with each other because the nonlinear dynamics can be described only within few theoretical models dodd . Thus there exists a possibility to predict novel effects in realistic physical systems via modeling similar processes in simple hypothetical systems, namely chains of coupled nonlinear oscillators could be used for this purpose. For instance, as recently has been discovered by Geniet and Leon leon  nonlinear supratransmission phenomenon takes place in discrete sine-Gordon lattice, this means that driving harmonically and continuously one end of the lattice with frequencies within a band gap, there is no energy flow through the lattice for low amplitude driving, while above definite driving amplitude threshold a sudden energy flow takes place. This nontrivial effect has been explained by means of the direct soliton creation at the end of the lattice, in other words the sudden energy flow occurs when the driving adjusts the internal oscillations of the localized object. It was also noted there the possibility of the existence of similar mechanism of gap soliton generation in photonic band gap materials. It should be especially mentioned that nonlinear supratransmission has been detected not only making numerical simulations for model system of discrete sine-Gordon lattice, but it is also experimentally realized on a mechanical pendulums chain driven at one end at band gap frequencies leon .

The present letter aims to analyze whether the similar scenario takes place in case of discrete nonlinear Schrödinger (DNLS) equation and then make the predictions concerning the corresponding nonlinear processes in coupled optical waveguide arrays christo . The experimental conditions are suggested for which optical waveguide array becomes transparent with respect to the beam injected into the single boundary waveguide if beam’s intensity exceeds certain threshold (see Fig 1). This effect is opposite to the ordinary case when for low intensities (linear regime) the light injected into the single waveguide spreads to other waveguides, while in nonlinear case the light is trapped into several neighboring waveguides leading thus to the spatial discrete optical breather creation (see e.g. most recent experimental papers on the subject erice1 ). It should be just mentioned that longitudinal space dimension in optical waveguide array plays a role of a time variable and one should especially care about this suggesting real experiments on waveguide arrays.

Let us start form the consideration of boundary driven DNLS equation which could be written in the following form (j=1, $`\mathrm{}`$ , N):

$$i\frac{\psi _j}{z}+\psi _{j+1}+\psi _{j1}+2|\psi _j|^2\psi _j=0;\psi _0=Ae^{i\mathrm{\Delta }z}.$$
(1)
Here $`z`$ variable stands for the time, $`\mathrm{\Delta }`$ and $`A`$ are driving frequency and amplitude of the boundary. Initial condition reads as $`\psi _j(0)=0`$ and it is supposed that the driving is turned on adiabatically, e.g. $`A=A_0\left[1\mathrm{exp}(z/\tau )\right]`$ in order to avoid the appearance of the perturbations from the initial shock. In simulations it is taken $`\tau =10`$ and nonlinear dynamics is monitored up to the time scales $`10^4`$, thus in stationary regime ($`z\tau `$) one has $`A=A_0`$. The damping also has been applied at the right end of the waveguide array in order to suppress edge reflection. Note that in the absence of driving the sum of intensities $`_j^N|\psi _j|^2`$ is a conserved quantity, thus the nonlinear dynamics across the array could be described via intensity flux $`J_j`$ through the site $`j`$:

$$J_j=i\left(\psi _j\psi _{j+1}^{}\psi _j^{}\psi _{j+1}\right)$$
(2)
where the intensity and intensity flux at site $`j`$ are connected with each other via discrete continuity condition $`d|\psi _j|^2/dz+\left(J_jJ_{j1}\right)=0`$.

The numerical simulations have been performed choosing different values of boundary driving parameters $`\mathrm{\Delta }`$ and $`A`$. From the numerical experiments it follows that boundary driving leads to the perturbation of all sites if driving frequency $`\mathrm{\Delta }`$ is located within the linear phonon band $`2<\mathrm{\Delta }<2`$, i.e. there is a nonzero intensity flux for any driving amplitudes (see Fig. 2). On the other hand if the driving frequency is in upper band gap, particularly $`\mathrm{\Delta }>2`$, for low driving amplitudes only several neighboring sites are excited and intensity flux to remote sites is zero. The energy starts to flow only if the driving amplitude exceeds certain threshold (see Figs. 3 and 4). Note that if the driving frequencies are within a lower ban gap there is no intensity flow for any driving amplitudes. Now it is time to analyze the mechanisms for that effect.

As far as a boundary driving is applied it is natural to expect that localized solutions will excite. This statement is in full accordance with the consideration of similar process in discrete sine-Gordon type models where the nonlinear supratransmission has been discovered leon . Thus one can assume that nonzero intensity flux will appear when boundary driving excites moving localized solutions. It is easy to derive a semi-discrete approximate envelope soliton solution substituting ansatz $`\psi _j=\mathrm{\Psi }(j)\mathrm{exp}\{i(\beta z\chi j)\}`$ into the DNLS equation (1). Then assuming that the envelope $`\mathrm{\Psi }(j)`$ varies smoothly along the lattice and expanding $`\mathrm{\Psi }(j\pm 1)`$ about the site $`j`$ we get the following approximate one soliton solution (see for the details of similar derivation e.g. in Ref. rani ):

$$\psi _j=\frac{|\psi _j|_{max}}{\mathrm{cosh}\left[|\psi _j|_{max}(jVz)\right]}e^{i(\beta z\chi j)}$$
(3)
with a following nonlinear dispersion relation for the carrier wave of the envelope soliton ($`\chi `$ varies from 0 to $`2\pi `$)

$$\beta =2\mathrm{cos}\chi +|\psi _j|_{max}^2,$$
(4)
and $`V=\beta /\chi =2\mathrm{sin}\chi `$ is a soliton’s group velocity. Note that the assumption that soliton envelope varies smoothly along the lattice puts the following restriction on the soliton amplitude $`|\psi _j|_{max}^21`$.

It is expected that intensity flux appears in the system if driving adjusts the nonlinear dispersion relation (4), i.e. $`\beta =\mathrm{\Delta }`$ and $`|\psi _j|_{max}=A_0`$. In other words, flux is nonzero only if one can find such $`\chi `$ that the following condition is fulfilled $`\mathrm{\Delta }=2\mathrm{cos}\chi +A_0^2`$. Therefore for inband driving $`2<\mathrm{\Delta }<2`$ the nonzero flux appears even for very low driving amplitudes $`A_0`$, while in upper band gap

$$\mathrm{\Delta }>2$$
(5)
there exists certain amplitude threshold

$$A_{th}=\sqrt{\mathrm{\Delta }2}$$
(6)
below which ($`A_0<A_{th}`$) there is no intensity flux into the system. Instead, only several sites are excited and that pattern could be described lepri  by a static breather solution of Eq. (1). This solution could be derived from the general one (3) for zero velocity $`V=2\mathrm{sin}\chi =0`$, i.e. $`\chi =0`$ and requiring the fulfillment of boundary condition $`\psi _0=A_0\mathrm{exp}(i\mathrm{\Delta }z)`$ in stationary state ($`z\tau `$):

$$\psi _j=\frac{\sqrt{\mathrm{\Delta }2}}{\mathrm{cosh}\left[(j+x_0)\sqrt{\mathrm{\Delta }2}\right]}e^{i\mathrm{\Delta }z},$$
(7)
where $`x_0=\text{acosh}\left[\sqrt{\mathrm{\Delta }2}/A_0\right]/\sqrt{\mathrm{\Delta }2}`$.

The results of numerical simulations are fully explained by the above consideration. In inset of Fig. 3 it is presented the comparison between numerical experiments on Eq. (1) and analytical formula (6) for the driving threshold above which ($`A_0>A_{th}`$) a nonzero intensity flux appears in the system. For driving frequencies $`\mathrm{\Delta }`$ close to 2 there is a perfect agreement, but this agreement becomes worse for larger driving frequencies. The point is that for driving frequencies sufficiently larger than 2 threshold amplitudes become comparable with unity according to the relation (6). But for such amplitudes the continuum envelope approximation (3) is invalid, moreover, as numerical simulations show large amplitude excitations are trapped by the lattice, (see Fig. 5), as a result localizations do not move and intensity flux becomes zero. Thus the band gap transmission effect exists if there are moving solutions in the system. As a result in case of DNLS equation the discovered phenomenon is observable for driving frequencies $`2<\mathrm{\Delta }<2.09`$.

Now let us discuss how the obtained results could be applied to describe nonlinear transmission processes in the system of coupled optical waveguides. For the purpose to realize a band gap driving it is suggested (see Fig. 1) to inject a beam into boundary waveguide with linear refractive index $`n_0`$ larger than refractive index $`n`$ of other waveguides forming the array. Let us introduce a linear coupling constant between array waveguides as $`Q`$, while the coupling between boundary ($`j=0`$) and first waveguides is defined as $`Q_0`$. Besides that, let us suppose for simplicity that the nonlinear refractive index (kerr nonlinearity) in $`j=0`$ waveguide is equal to zero and onsite nonlinear refractive index in array waveguides is $`D`$. Thus the wave envelopes in waveguides could be described by a set of the following equations ($`j=2,\mathrm{},N`$):

$`i{\displaystyle \frac{\psi _0}{z}}+{\displaystyle \frac{\omega }{c}}n_0\psi _0+Q_0\psi _1=0`$
$`i{\displaystyle \frac{\psi _1}{z}}+{\displaystyle \frac{\omega }{c}}n\psi _1+Q_0\psi _0+Q\psi _2+D|\psi _1|^2\psi _1=0`$ (8)
$`i{\displaystyle \frac{\psi _j}{z}}+{\displaystyle \frac{\omega }{c}}n\psi _j+Q(\psi _{j+1}+\psi _{j1})+D|\psi _j|^2\psi _j=0.`$
where $`\omega `$ is a carrier wave frequency and $`c`$ is a light velocity. Last equation from the set (8) is a well known representation of infinite waveguide array by the DNLS equation christo  while first two equations describe influence of boundary waveguide (with different linear refraction index) on the semi infinite array (see also Ref. kivshar ). After the appropriate rescaling

$`\psi _j=\psi _j^{}e^{iz\omega n/c}\sqrt{2Q/D}\text{for}j=1,\mathrm{},N;`$ (9)
$`\psi _0=\psi _0^{}e^{iz\omega n/c}\left(Q/Q_0\right)\sqrt{2Q/D};z=z^{}/Q`$
Eqs. (8) obtain simpler form ($`\mathrm{\Delta }\omega \left(n_0n\right)/Qc`$):

$`i{\displaystyle \frac{\psi _0^{}}{z^{}}}+\mathrm{\Delta }\psi _0^{}+{\displaystyle \frac{Q_0^2}{Q^2}}\psi _1^{}=0`$ (10)
$`i{\displaystyle \frac{\psi _j^{}}{z^{}}}+(\psi _{j+1}^{}+\psi _{j1}^{})+|\psi _j^{}|^2\psi _j^{}=0.`$
($`j=1,\mathrm{},N`$) which reduces to the boundary driven DNLS (1) with boundary condition $`\psi _0^{}=\psi _0^{}(0)\mathrm{exp}(i\mathrm{\Delta }z^{})`$ in the limit $`(Q_0/Q)0`$ and therefore for $`Q_0/Q1`$ one can use the results derived for the case of boundary driven DNLS (5) and (6). Particularly, for

$$\mathrm{\Delta }\omega \left(n_0n\right)/Qc>2$$
(11)
the localized excitations (3) or (7) form with propagation constant located in the upper band gap $`\beta =\mathrm{\Delta }>2`$. Thus if the injected intensity in the boundary is below the threshold

$$|\psi _0(0)|_{th}^2=\frac{2Q^3}{DQ_0^2}|\psi _0^{}(0)|_{th}^2\frac{2Q^2}{DQ_0^2}\left[\frac{(n_0n)\omega }{c}2Q\right]$$
(12)
one has a static breather solution (7) and intensity flux to the array waveguides is equal to zero, while above the threshold energy transmission begins via gap solitons (3).

The above expression for the threshold (12) becomes exact in the limit $`(Q_0/Q)0`$. For small but nonzero $`(Q_0/Q)`$ the beam intensity in the boundary waveguide could be considered as almost constant quantity irrespective to the spread of energy in nonlinear band gap transmission regime, because according to the rescaling beam intensity in the boundary waveguide is $`Q/Q_0`$ times larger than the amplitude of propagating soliton through the array. In case when $`Q_0`$ becomes comparable with $`Q`$ almost all intensity in boundary waveguide is needed to form a soliton. As a result the intensity in boundary waveguide sharply decreases and therefore much larger threshold intensity is needed \[than given by relation (12)\] to develop the nonlinear transmission in band gap regime. In Fig. 6 it is presented the picture for the case $`Q_0/Q=0.1`$. As seen the beam intensity in boundary waveguide is above the threshold at the origin $`z^{}=0`$ and it produces gap solitons causing the band gap transmission. But this process itself causes the decrease of beam intensity in the boundary waveguide, the intensity after the creation of several solitons goes below the threshold and transmission process is not observable for large $`z^{}`$.

As it has been mentioned above nonlinear band gap transmission regime holds if one has low amplitude solitons. Large amplitude solitons tend to pin and energy transfer becomes much less effective. However, this is true only in case of DNLS equation where only onsite nonlinearities are taken into account. As shown recently erice2  considering also the terms describing also intersite nonlinearities one has moving breather solutions even at large excitations amplitudes. High intensity moving breathers have been also detected on the recent experiments erice1 . Therefore the numerical simulations have been undertaken adding to DNLS equation also terms with intersite nonlinearities. In this case the energy transfer via the moving breathers take place even for large excitation amplitudes and as a result optical transparency regime is observable in whole range of $`\mathrm{\Delta }>2`$.

Summarizing it should be noted again that the novel scenario of nonlinear band gap transmission in optical waveguide arrays is predicted and simple experimental setup is suggested for its realization. Suggested experimental setup would serve also for generation of optical gap solitons propagating across the waveguide array.

Finally I would like to thank all participants of NATO advanced research workshop at Erice (Italy, 2003) for inspiring discussions. Special thank go to Yuri Kivshar, Stefano Lepri, Stefano Ruffo and Lasha Tkeshelashvili for suggestions. The work is supported by USA CRDF award No GP2-2311-TB-02.


text:
# Phase separation in the bosonic Hubbard model with ring exchange

## Abstract

We show that soft core bosons in two dimensions with a ring exchange term exhibit a tendency for phase separation. This observation suggests that the thermodynamic stability of normal bose liquid phases driven by ring exchange should be carefully examined.

Interest in ring exchange interactions in quantum many-body systems has a long history, both theoretically and experimentallyThouless65 . Recently, the ring exchange interaction has been invoked in an effort to understand various aspects of high temperature superconductivity. While the Heisenberg model alone provides a rather accurate picture of magnetic excitations in the parent compounds of the cuprate superconductorsManousakis91 , estimates of the magnitude of the ring exchange term are as high as one quarter of the exchange couplingColdea01 ; Katanin02 ; MullerHartmann02  and it therefore has been of interest to understand how this term might modify magnetic propertiesRoger89 ; Honda93 ; Lorenzana99 ; Coldea01 ; Matsuda00 . Ring exchange interactions have also been suggested as a likely candidate to reconcile the properties of the underdoped pseudogap regime. The basic picture is that the ring exchange interaction can give rise to a new normal “Bose metal” phase at zero temperature in which there are no broken symmetries associated with superfluidity or charge density wave phases, and in which the compressibility is also finiteFisher02 .

With these motivations partly in mind, Sandvik et alSandvik02  studied the phase diagram of the two-dimensional spin-$`1/2`$ $`XY`$ model with spin exchange interaction on a square lattice,

$$H=J\underset{\mathrm{𝐢𝐣}}{}B_{\mathrm{𝐢𝐣}}K\underset{\mathrm{𝐢𝐣𝐤𝐥}}{}P_{\mathrm{𝐢𝐣𝐤𝐥}}$$
(1)
where

$`B_{\mathrm{𝐢𝐣}}`$ $`=`$ $`S_𝐢^+S_𝐣^{}+S_𝐢^{}S_𝐣^+=2(S_𝐢^xS_𝐣^x+S_𝐢^yS_𝐣^y),`$
$`P_{\mathrm{𝐢𝐣𝐤𝐥}}`$ $`=`$ $`S_𝐢^+S_𝐣^{}S_𝐤^+S_𝐥^{}+S_𝐢^{}S_𝐣^+S_𝐤^{}S_𝐥^+`$ (2)
and $`\mathrm{𝐢𝐣}`$ denotes nearest neighbours and $`\mathrm{𝐢𝐣𝐤𝐥}`$ are sites at the corners of a plaquette. As is well known, for $`K=0`$ this model is exactly equivalent to the hard core bosonic Hubbard model, at half filling, with no interactions apart from the constraint on site occupations, and it has only a superfluid phase. Sandvik et alSandvik02  studied the phase diagram as a function of $`K`$ using the Stochastic Series Expansion algorithmSandvik91 . They found that as $`K`$ increases, the superfluid density, $`\rho _s`$, decreases up to a critical value, $`K_c`$, where a phase transition takes place, and $`\rho _s`$ goes to zero with long range order appearing in the momentum $`(\pi ,0)`$ and $`(0,\pi )`$ channels of the plaquette-plaquette correlation function. This indicates the existence of a phase in which plaquettes with large and small values of ring exchange alternate with a striped pattern across the lattice. This phase is also an incompressible insulator. For yet larger $`K`$, charge density wave order is established in which the site occupations are alternatingly large and small.

It is believed Fisher02  that when the hardcore constraint is relaxed, the striped plaquette phase might evolve into a normal compressible conducting “Bose metal” in which none of the order parameters mentioned above is non-zero. This suggestion leads us to study here the phase diagram of the soft core bosonic Hubbard model at half filling with ring exchange interaction,

$`H`$ $`=`$ $`t{\displaystyle \underset{𝐢,𝐣}{}}(a_𝐢^{}a_𝐣+a_𝐣^{}a_𝐢)+U{\displaystyle \underset{𝐢}{}}n_𝐢(n_𝐢1)`$
$`+K{\displaystyle \underset{\mathrm{𝐢𝐣𝐤𝐥}}{}}(a_𝐢^{}a_𝐣a_𝐤^{}a_𝐥+a_𝐢a_𝐣^{}a_𝐤a_𝐥^{})`$
where the destruction and creation operators satisfy $`[a_𝐢,a_𝐣^{}]=\delta _{\mathrm{𝐢𝐣}}`$, $`n_𝐢=a_𝐢^{}a_𝐢`$ is the number operator at site $`𝐢`$ and $`U`$ is the onsite interaction strength. For our quantum Monte Carlo (QMC) simulations we used the World Line algorithm with four-site decouplingLoh85 . We verified our code for $`K=0`$ by comparing with existing results for hard and soft core bosons with and without near and next near neighbor interactions. For the $`K0`$ case, we compared with the hard core results of Sandvik02 .

Before discussing results for the full many-body system, it is interesting to study the behaviour of two bosons, since the formation of a bound state is closely related to the issue of phase separation. Ring exchange, like an attractive potential, favors proximity of two bosons, since the action of such a term is nonzero only when the two bosons live on the same plaquette. In Fig. 1 we show the average separation $`\mathrm{\Phi }_0|r^2|\mathrm{\Phi }_0`$ of two bosons, normalized to the number of sites $`L^2`$ on an $`L`$x$`L`$ lattice. As $`K`$ increases, there is a crossover at $`K/t3`$ from a regime where the boson separation grows linearly with system size, so that the normalized separation is size independent, to one in which the boson separation does not grow with system size.

While we do not show the associated data, plots of $`\mathrm{\Phi }_0|r^2|\mathrm{\Phi }_0`$ for different soft-core repulsion $`U`$ reveal that the average separation is insensitive to the value of $`U`$. This will obviously be true in the unbound regime at small $`K`$, since the density is so dilute. It is less clear that this should be so in the bound regime at large $`K`$. However, as can be seen from the data in Fig. 1, the radius of the bound state is several lattice spacings ($`r^20.1L^2`$ whence $`r0.3L`$), so here too the effect of $`U`$ is expected to be relatively small.

Fig. 1 suggests that there might be a tendency for ring exchange to cause the bosons to clump together, and, in an extreme scenario, to undergo phase separation. However, at densities higher than the dilute two boson case, this effect is opposed by the repulsion $`U`$. The focus of this paper is to examine this competition and determine the phase diagram of the soft core case as a function of both $`U`$ and $`K`$ at half filling.

The most straightforward indication of phase separation comes from a real space image of the boson density during the course of a simulation. Fig. 2 shows the average density distributionnote2  for $`L=16`$, $`U=4`$ and $`K=2.5`$. We see indications that the bosons undergo phase separation: At less than quarter filling (top panel) the bosons clump together into a compact region of high density. At densities above quarter filling, on a lattice with periodic boundary conditions, the number of occupied plaquettes is largest (and hence the ring exchange energy is most negative) for a configuration where the bosons stretch out in stripe across the lattice (bottom panel).

We will now demonstrate that phase separation is characteristic of a large portion of the $`KU`$ phase diagram by examining the density-density correlation function and its associated structure factor, fixing $`U`$ (or $`K`$) and scanning $`K`$ (or $`U`$). As Fig. 2 illustrates, if the bosons phase separate, they may form a structure in which a set of contiguous sites of about half the system size will have appreciable boson occupation. The other half of the lattice is essentially empty. Therefore, if one examines the structure factor of the density-density correlation function,

$$S(k_x,k_y)\frac{1}{L^2}\underset{𝐫}{}C(𝐫)\mathrm{e}^{i𝐫.𝐤}$$
(4)
with

$$C(𝐫)=\frac{1}{L^2}\underset{𝐫^{}}{}n(𝐫^{})n(𝐫+𝐫^{}),$$
(5)
one should observe a peak in $`S`$ at small momentum, e.g. $`(2\pi /L,0)`$, $`(0,2\pi /L)`$ or $`(2\pi /L,2\pi /L)`$ depending on the precise orientation of the clumpnote3 . By looking at the sum of the density structure factor at these three smallest momentum values we are sensitive to phase separation regardless of whether it occurs in a puddle of roughly circular shape (Fig. 2a) or in some more elongated pattern (Fig. 2b). Fig. 3 clearly shows this behaviour: For $`K<2`$ at $`U=4`$, $`S`$ is very small at the relevant momenta. For $`K>2`$, phase separation sets in. Data for $`16\times 16`$ and $`24\times 24`$ lattices are shown and their agreement indicates that this phase separation is not a finite lattice effect. The critical value of $`K`$ grows roughly linearly with $`U`$.

It is also interesting to understand the behaviour of the superfluid density $`\rho _s`$. One does not necessarily expect $`\rho _s`$ to vanish when phase separation occurs. In fact, as is well known, $`\rho _s0`$ for the soft-core boson Hubbard model at all fillings, including commensurate density, if $`U/t`$ is sufficiently small. Similarly, here it is possible that $`\rho _s`$ can survive in the dense region of the phase separated latticeFisher89 ; Batrouni90 . Our simulations show that when phase separation first occurs, the populated region forms a band that spans the whole system. The bosons may then delocalize along that band, maintaining an (anisotropic) superfluid density. We have found that when the bosons form such a band, the plaquette-plaquette structure factor is also anisotropic and has long range correlations along the direction of the band. As $`K`$ is increased further, the populated region of the lattice takes the form of an island. In such a case, the system may not be considered a superfluid in that one cannot establish superflow across the system. However, the bosons may still be delocalized over the extent of the islandnote4 . Fig. 4 shows $`\rho _s`$ versus $`K`$.

By making scans like those of Fig. 3 at several values of $`U`$, we construct the phase diagram which we show in Fig. 5. Above the solid line, the system undergoes phase separation.

How do these results connect with the previous studies in the hard-core limit? There we know that a phase transition occurs at $`K_c/J7.9`$ from a superfluid to a striped plaquette phaseSandvik02 . We show in Fig. 6 the behaviour of the boson density $`\rho `$ as a function of chemical potential $`\mu `$ in this hard-core limit. The jump in $`\mu `$ across half-filling $`\rho =0.5`$ shows that the plaquette ordered phase has a gap to the addition of bosons (“charge excitations”), a result which is in agreement with Sandvik etal Sandvik02 . The slope of the $`\rho `$ versus $`\mu `$ is the compressibility $`\kappa `$. Consequently, if this curve “bends backwards”, the system is thermodynamically unstable and undergoes phase separationBatrouni00 . While the data are not conclusive, we do see hints of an instability. For $`|\rho 0.5|>0.008`$, the slope is finite and corresponds to a normal superfluid. For $`|\rho 0.5|<0.008`$, the slope is either very large, or perhaps negative, indicative of phase separation.

So far we have addressed mostly the half filled case, and lower densities. It is of course of interest to examine higher fillings where the effect of $`U`$ will be expected to discourage phase separation. We have done simulations for $`\rho =1`$ and $`\rho =1.5`$ and found in both cases that despite the higher densities, phase separation still sets in above a critical value $`K>K_c`$ of the ring exchange energy scale. This result is not so surprising since the on-site repulsion $`U`$ and the ring exchange term both scale with density as $`\rho ^2`$. (Each of the four ring exchange creation/destruction operators picks up a factor of $`\sqrt{n}`$ when acting on a site with occupation $`n`$.)

In conclusion, we have shown that a sufficiently large ring exchange energy can lead to a thermodynamic instability and phase separation. We determined the critical $`K`$ as a function of the soft core repulsion $`U`$ for a half-filled lattice and found roughly $`K_cU/2`$. We conclude that the soft core boson Hubbard model does not exhibit a normal Bose metal phase. The bosons are either superfluid or undergo phase separation. This phase separation also takes place when the hopping parameter vanishes, $`t=0`$, a limit examined for the quantum phase model in referenceFisher02  but which did not find phase separation.

Finally, let us comment on the implications of our work for the phase diagram of the spin-1/2 quantum Heisenberg model with a ring exchange term. The kinetic energy term in the hard-core boson Hubbard model maps onto $`J(S_𝐢^xS_𝐣^x+S_𝐢^yS_𝐣^y)`$ with exchange constant $`J=2t`$. At the value $`U=4t`$ in our soft core model, double occupancy is already very rare at half-filling, and hence we are almost in the hard-core limit. The value of the ring exchange energy scale required to drive phase separation for this $`U`$ is $`K2t`$, or in other words, $`KJ`$. To replicate the near-neighbor coupling of the $`z`$ components of spin in the Heisenberg model we must include a near-neighbor repulsion in the bose-Hubbard model, a term which clearly would suppress phase separation. Thus we expect ring exchange to have the potential to drive phase separation in the Heisenberg model only for $`K`$ considerably greater than $`J`$.

Note added: After the completion of this work a preprint by Melko, Sandvik and Scalapino (cond-mat/0311080) appeared where the phase diagram of the doped hardcore system with ring exchange was determined.

Acknowledgments: We acknowledge useful conversations with T. C. Newman. This work was supported by NSF-CNRS cooperative grant #12929, NSF-DMR-0312261 and NSF-INT-0124863.