Abstract

We present two-dimensional (2D) hydrodynamic simulation of rotating galactic winds driven by radiation. We study the structure and dynamics of the cool and/or warm component ( T ≃ 10 4 K) which is mixed with dust. We have taken into account the total gravity of a galactic system that consists of a disc, a bulge and a dark matter halo. We find that the combined effect of gravity and radiation pressure from a realistic disc drives the gas away to a distance of ∼5 kpc in ∼37 Myr for typical galactic parameters. The outflow speed increases rapidly with the disc Eddington parameter Γ 0 (=κ I /(2 cG Σ)) for Γ 0 ≥ 1.5. We find that the rotation speed of the outflowing gas is ≲100 km s −1 . The wind is confined in a cone that mostly consists of low angular momentum gas lifted from the central region.

1 Introduction

Many galaxies are observed to have moving extraplanar gas, generally termed as galactic superwinds (see Veilleux, Cecil & Bland-Hawthorn 2005 for a recent review). Initial observations showed the Hα emitting gas above the plane of M82 (e.g. Lynds & Sandage 1963 ). The advent of X-ray astronomy established yet another phase of galactic outflows, namely the hot plasma, emitting X-rays in the temperature range 0.3–2 keV ( Strickland et al. 2004 ). Also recent observations have revealed the existence of molecular gas in these outflows ( Walter, Weiss & Scoville 2002 ; Veilleux, Rupke & Swaters 2009 ). Earlier observations were limited to local dwarf starburst galaxies that showed these winds. However, in recent years, the observations of outflows in Ultra-Luminous Infra-red Galaxies (ULIGs) have extended the range of galaxies in which outflows are found ( Rupke, Veilleux & Sanders 2002 , 2005 ; Martin 2005 ).

On the theoretical side, there have been speculations on winds from starburst galaxies ( Burke 1968 ; Johnson & Axford 1971 ; Mathews & Baker 1971 ). In these models the large-scale winds are a consequence of energy injection by multiple supernovae (SNe; Larson 1974 ; Chevalier & Clegg 1985 ; Dekel & Silk 1986 ; Heckman 2002 ). In the context of the multiphase structure of the outflows, the results of these theoretical models are more relevant for the X-ray emitting hot wind. On the other hand, observations of the cold outflows are better explained by the radiation driving ( Martin 2005 ; Murray, Quataert & Thompson 2005 ).

If only the Thompson scattering is considered, then radiation from galaxies does not seem to be a reasonable wind-driving candidate because opacities would be small; however one should consider that these winds are heavily enriched. Murray et al. (2005) proposed a wind-driving mechanism based on the scattering of dust grains by the photons from the galaxy (see also Chiao & Wickramasinghe 1972 ; Davies et al. 1998 ). This mechanism can be quite effective since the opacities in dust–photon scattering can be of the order of one hundred cm 2 g −1 and gas in turn, being coupled with the dust, is driven out of the galaxy if the galaxy possesses a certain critical luminosity. Bianchi & Ferrara (2005) argued that dust grains ejected from galaxies by radiation pressure can enrich the intergalactic medium (IGM). Nath & Silk (2009) then described a model of outflows with radiation and thermal pressure in the context of outflows from Lyman break galaxies observed by Shapley et al. (2003) . Murray et al. (2011) have also described a similar model in which radiation pressure is important for the first few million years of the starburst phase, after which SN-heated hot gas pushes the outflowing material. Sharma & Nath (2011) have also shown that radiation pressure is important for outflows from high-mass galaxies with a large star formation rate (SFR) (with vc ≥ 200 km s −1 , SFR ≥ 100 M yr −1 ), particularly in ULIGs.

In this paper, we study the effect of radiation pressure in driving cold and/or warm gas outflows from disc galaxies with numerical simulations. Recently, Sharma, Nath & Shchekinov (2011) calculated the terminal speed of such a flow along the pole of a disc galaxy, taking into account the gravity of the disc, the stellar bulge and the dark matter halo. They determined the minimum luminosity (or, equivalently, the maximum mass-to-light ratio of the disc) to drive a wind, and also showed that the terminal speed lies in the range of 2–4 vc (where vc is the rotation speed of the disc galaxy), consistent with observations ( Martin 2005 ; Rupke et al. 2005 ), and the ansatz used by numerical simulations in order to explain the metal enrichment of the IGM ( Oppenheimer & Davé 2006 ). We investigate further the physical processes for a radiation-driven wind. Rotation is yet another aspect of the winds that we address in our simulation. As the wind material is lifted from a rotating disc, it should be rotating inherently as seen in observations as well ( Sofue et al. 1992 ; Seaquist & Clark 2001 ; Walter et al. 2002 ; Greve 2004 ; Wesmoquette et al. 2009 ).

Previous simulations of galactic outflows have considered the driving force of a hot interstellar medium (ISM) energized by the effects of SNe ( Kohji & Ikeuchi 1988 ; Tomisaka & Bregman 1993 ; Suchkov et al. 1994 , 1996 ; Mac Low & Ferrara 1999 ; Strickland & Stevens 2000 ; Fragile, Murray & Lin 2004 ; Cooper et al. 2008 ; Fujita et al. 2009 ). However, the detailed physics of a radiatively driven galactic outflow is yet to be studied with a simulation. In this work, we study the dynamics of an irradiated gas above an axisymmetric disc galaxy by using hydrodynamical simulation. Recently Hopkins, Quataert & Murray (2012) have explored the relative roles of radiation and SNe heating in galactic outflows, and studied the feedback on the star formation history of the galaxy. Our goal here is different in the sense that we focus on the structure and dynamics, particularly the effect of rotation, of the wind. In order to disentangle the effects of various processes involved, we intentionally keep the physical model simple. For example, we begin with a constant density and surface brightness disc, then study the effect of a radial density and radiation profile, and finally introduce rotation of the disc, in order to understand the effect of each detail separately, instead of performing one single simulation with many details put together.

The main driving force is radiation force and the containing force is due to gravity. We take the system to be composed of three components: a disc, a bulge and a dark matter halo. We describe the forces due to these three constituents below. We take a thin galactic disc and a spherical bulge. All these forces are given in cylindrical coordinates because we solve the fluid equations in cylindrical geometry.

2.1 Gravitational field from the disc

Consider a thin axisymmetric disc in r φ plane with surface mass density Σ( r ). As derived in Appendix A, the vertical and radial components of gravity due to the disc material at a point Q above the disc with coordinates ( r , 0, z ) are given by

1
The azimuthal coordinate of Q is taken to be zero because of axisymmetry. The integration limit for φ′= 0 to .

We consider two types of disc in our simulations, one with uniform surface mass density and radius rd [uniform disc (UD)], and another with an exponential distribution of surface mass density [exponential disc (ED)] with a scale radius rs . The surface mass density of uniform surface density disc (i.e. UD) is

2
and in the case of a disc with exponentially falling density distribution (ED)
3
In case of UD ( equation 2 ), the integration limit would be r ′= 0 to rd , while for ED ( equation 3 ), the limits of the integration run from r ′= 0 to ∞. Numerically, this means we integrate up to a large value of r ′, increasing which will not change the gravitational field by any significant amount. We have chosen the Σs in such a way that the total disc mass remains the same for the UD or ED. Therefore,
4
In Fig. 1 , we plot the contours of gravitational field strength and its direction vectors due to a UD (left-hand panel), and that for the ED (right-hand panel). Interestingly, discs with the same mass but different surface density distributions produce different gravitational fields. For the UD the gravitational field is not spherical and the gravitational acceleration is maximum at the edge of the disc. On the other hand, the field due to ED is closer to spherical configuration with the maximum being closer to the centre of the disc and falling off outwards.

Figure 1

Magnitude of gravitational force of the (a) uniform disc (UD) and (b) exponential disc (ED) in colours with direction in arrows. Values are in the units of G Σ 0 (= 4.5 × 10 −9 cm s −2 ).

Figure 1

Magnitude of gravitational force of the (a) uniform disc (UD) and (b) exponential disc (ED) in colours with direction in arrows. Values are in the units of G Σ 0 (= 4.5 × 10 −9 cm s −2 ).

2.2 Bulge and the dark matter halo

We consider a bulge with a spherical mass distribution and constant density, with mass Mb and radius rb . The radiation force due to the bulge is negligible as it mostly hosts the old stars. The gravitational force of the bulge is given by

5

6
where .

We consider an NFW halo ( Navarro, Frenk & White 1997 ) with a scaling with disc mass as given by Mo, Mao & White (1998 , hereafter MMW98) where the total halo mass is ∼20 times the total disc mass. The mass of an NFW halo has the following functional dependence on R :

7
where . Here, ρ crit is the critical density of the Universe at present epoch, Rs is the scale radius of the NFW halo and R200 is the limiting radius of virialized halo within which the average density is 200ρ crit . This mass distribution corresponds to the following potential:
8
The gravitational force due to the dark matter halo is therefore given by
9
The net gravitational acceleration is therefore given by
10
The gravitational field for both bulge and halo is spherical in nature, although the field due to the bulge maximizes at rb . However, the net gravitational field will depend on the relative strength of the three components. In Fig. 2 (left-hand panel), we plot the contours of total gravitational field strength due to the bulge, the halo and a UD. The non-spherical nature of the gravitational field is evident. A more interesting feature appears due to the bulge gravity. The net gravitational intensity maximizes in a spherical shell of radius rb (= 0.2 Lref ; see Section 3.1 ). Therefore, there is a possibility of piling up of outflowing matter at around a height zrb near the axis. In the right-hand panel of Fig. 2 , we present the contours of net gravitational field due to an embedded ED within a halo and a bulge.

Figure 2

Total gravitational force of the (a) UD and (b) ED in colours with direction in arrows. The values are in the same units as in Fig. 1 .

Figure 2

Total gravitational force of the (a) UD and (b) ED in colours with direction in arrows. The values are in the same units as in Fig. 1 .

2.3 Radiation from disc and the Eddington factor

We treat the force due to radiation pressure as it interacts with charged dust particles that are assumed to be strongly coupled to gas by Coulomb interactions and which drags the gas with it. The strength of the interaction is parametrized by the dust opacity κ which has the units cm 2 gm −1 .

Gravitational pull on the field point Q ( R , Z ) due to the disc point P ( r ′, φ′, 0) is along the direction (see Appendix A). The difference in computing the radiation force arises due to the fact that one needs to account for the projection of the intensity at Q (for radiation force from more complicated disc, see Chattopadhyay 2005 ). For a disc with surface brightness I ( r ), we can find the radiation force by replacing G Σ( r ′) in equation (1) by I ( r ′)κ/ c , and taking into account the projection factor . Similar to the disc gravity, the net radiation force Frad at any point will have the radial component ( F rad, r ) and the axial component ( F rad, z ) and are given by

11

12
Since we have two models for disc gravity, we also consider two forms of disc surface brightness.
13
and
14
If the two disc types are to be compared for identical luminosity, then one finds
15
The disc Eddington factor is defined as the ratio of the radiation force and the gravitational force (MQT05). In spherical geometry this factor is generally constant at each point because both gravity and radiation have an inverse square dependence on distance. Although in the case of a disc, the two forces have different behaviour, we can still define an Eddington parameter as . In this case this parameter depends on the coordinates r , φ, z of the position under consideration. We can, however, define a parameter whose value is the Eddington factor at the centre of the disc, i.e.
16
If Γ 0 = 1, then the radiation and gravity of the disc will cancel each other at the centre of the disc. We will parametrize our results in terms of Γ 0 . Therefore, the components of the net external force due to gravity and radiation is given by
17
In Fig. 3 , we plot the contours of radiative acceleration from a UD, and the same from an ED. There is a significant difference between the radiation field above an ED and that above a UD. The radiation field from a UD is largely vertical for small radii, but starts to diverge at the disc edge, at rrd . One can therefore expect that for high enough I , the wind trajectory will diverge. In the case of ED, the radiation field above the inner portion of the disc is strong and decreases rapidly towards the outer disc.

Figure 3

Magnitude of force due to radiation from the (a) UD and (b) ED for Γ 0 = 0.5, with arrows for direction.

Figure 3

Magnitude of force due to radiation from the (a) UD and (b) ED for Γ 0 = 0.5, with arrows for direction.

3 Numerical Method

The hydrodynamic equations have been solved in this paper by using the Total Variation Diminishing (TVD) code, which has been quite exhaustively used in cosmological and accretion disc simulations (see Ryu et al. 1993 , 1995 ; Kang et al. 1994 ; Molteni, Ryu & Chakrabarti 1996 ) and is based on a scheme originally developed by Harten (1983) . We have solved the equations in cylindrical geometry in view of the axial symmetry of the problem. This code is based on an explicit, second-order accurate scheme, and is obtained by first modifying the flux function and then applying a non-oscillatory first-order accurate scheme to obtain a resulting second-order accuracy (see Harten 1983 and Ryu et al. 1993 for details).

The equations of motion which are being solved numerically in the non-dimensional form is given by

18
where the state vector is
19
and the fluxes are
20
and the source function is given by
21

3.1 Initial and boundary conditions

We do not include the disc in our simulations and only consider the effect of disc radiation and total gravity on the gas being injected from the disc. We choose the disc mass to be Md = 10 11 M and assume it to be the unit of mass (i.e. Mref ). The units of length (i.e. Lref ) and velocity (i.e. vref ) are rd = 10 kpc and vc = 200 km s −1 , respectively. Therefore, the unit of time is tref = 48.8 Myr. We introduce a normalization parameter ξ such that , which turns out to be ξ= 1.08. Hence the unit of density is ρ ref = 6.77 × 10 −24 g cm −3 (∼4 mp cm −3 ). All the flow variables have been made non-dimensional by the choice of unit system mentioned above.

It is important to choose an appropriate initial condition to study the relevant physical phenomenon. We note that previous simulations of galactic outflows have considered a variety of gravitational potential and initial ISM configurations. For example, Cooper et al. (2008) considered the potential of a spherical stellar bulge and an analytical expression for disc potential, but no dark matter halo, and an ISM that is stratified in z -direction with an effective sound speed that is approximately five times the normal gas sound speed. Suchkov et al. (1994) considered the potential of a spherical bulge and a dark matter halo and an initial ISM that is spherically stratified. Fragile et al. (2004) considered a spherical halo and a z -stratified ISM. However, in a recent simulation of outflows driven by SNe from disc galaxies, Dubois & Teyssier (2008) found that the outflowing gas has to contend with infalling material from halo, which inhibits the outflow for a few Gyr. Fujita et al. (2004) also studied outflows from pre-formed disc galaxies in the presence of a cosmological infall of matter.

We choose a z -stratified gas to fill the simulation box, with a scale height of 100 pc. For the M 2 and M 3 case (of ED), we also assume a radial profile for the initial gas, with a scale length of 5 kpc. For the M 3 case, we further assume this gas to rotate with vφ decreasing with a scale height of 5 kpc. These values are consistent with the observations of Dickey & Lockman (1990) and Savage, Sembach & Lu (1997) for the warm neutral gas ( T ∼ 10 4 K) in Milky Way. We note that although the scale height for the warm neutral gas in our Galaxy is ∼400 pc at the solar vicinity, this is expected to be smaller in the central region because of the strong gravity due to bulge. The density of the gas just above the disc is assumed to be 0.1 particles cc −1 (0.025 in simulation units).

Furthermore, the adiabatic index of the gas is 5/3 and the gas is assumed initially to be at the same temperature corresponding to an initial sound speed cs (ini)=0.1 vref , a value that is consistent with the values in our Galaxy for the warm ionized gas with sound speed ∼18 km s −1 .

Our computation domain is rd × rd in the rz plane, with a resolution 512 × 512 cells. The size of individual computational cell is ∼20 pc. We have imposed a reflective boundary condition around the axis and zero rotational velocity on the axis. Continuous boundary conditions are imposed at r = rd and z = rd . The lower boundary is slightly above the galactic disc with an offset z0 = 0.01. We impose fixed boundary condition at lower z boundary. The velocity of the injected matter is vz ( r , z0 ) = v0 = 10 −5vref , and its density is given by

22
The density of the injected matter at the base (corresponding to 0.1 protons per cc).

For the case of ED with rotation (M 3 ), we assume for the injected matter to have an angular momentum corresponding to an equilibrium rotation profile. We show in Fig. 4 the rotation curves at z = 0 for all components (disc, bulge and halo) separately and the total rotation curve. We use the following approximation (shown by the thick red line in Fig. 4 ), which matches the total rotation curve,

23
We assume a bulge of mass Mb = 0.1 Mref and radius rb = 0.2 Lref . The scale radius for the NFW halo ( Rs ) is determined for a halo mass Mh = 20 Md , as prescribed by MMW98. The corresponding disc scale length is found to be rs ∼ 5.8 kpc, again using MMW98 prescriptions. Therefore we set the disc scale length for the ED case to be rs ∼ 0.58 Lref .

Figure 4

Rotation curves corresponding to the gravitational fields of an ED, bulge and halo are shown here in the units of vref [ = 200 km s −1 ], along with the total rotation curve. The approximation used in our simulation is shown by the thick red line.

Figure 4

Rotation curves corresponding to the gravitational fields of an ED, bulge and halo are shown here in the units of vref [ = 200 km s −1 ], along with the total rotation curve. The approximation used in our simulation is shown by the thick red line.

The above initial conditions have been chosen to satisfy the following requirements in order to sustain a radiatively driven wind as simulated here.

• The strong coupling between dust grains and gas particles requires that there are of the order of ∼ md / mp number of collisions between protons and dust grains of mass md ∼ 10 −14 g, for size m with density ∼3 g cm −3 . To ensure a sufficient number of collisions, the number density of gas particles should be cm −3 for Lref = 10 kpc.

• The time-scale for radiative cooling of the gas, assumed to be at T ∼ 10 4 K, is , where Λ∼ 10 −23 erg cm 3 s −1 ( Sutherland & Dopita 1993 ; table 6) for solar metallicity. The typical density filling up the wind cone in the realistic case (M 3 ) is ∼10 −3 to10 −4 cm −3 , which gives tcool ∼ 8–80 Myr and the dynamical time-scale of the wind is tref ∼ 50 Myr. Hence radiative cooling is marginally important, and we will address the issue of radiative cooling in a future paper.

• Radiative transfer effects are negligible since the total opacity along a vertical column of length Lref is κ( nmp ) Lref ∼ 0.003 for n ∼ 10 −3 cm −3 and κ∼ 100 cm 2 g −1 .

• The mediation of the radiation force by dust grains also implies that the gas cannot be too hot for the dust grains to be sputtered. The sputtering radius of grains embedded even in a hot gas of temperature T ∼ 10 5 K is m in a time-scale of 100 Myr ( Tielens et al. 1994 ), and this effect is not important for the temperature and density considered here.

3.2 Simulation set-up

We present three models with parameters listed in Table 1 . The initial conditions for all the models are described in Section 3.1 . The boundary condition is essentially the same, except that the mass flux into the computational domain from the lower z boundary depends on the type of disc. As has been mentioned in Section 3.1 , we keep the velocity of injected matter very low, vz ( r , z0 )= vz (ini)=10 −5vref , so that it does not affect the dynamics. The three models have been constructed by a combination of different values of three parameters Γ 0 , vφ and the distribution of the density in the disc. Model M 3 has been run for different values of Γ 0 to ascertain the effect of radiation.

Table 1

Models.

 Model name Γ 0 vφ Disc type M 1 2.0 0.0 UD M 2 2.0 0.0 ED M 3 2.0 1.0 ED
 Model name Γ 0 vφ Disc type M 1 2.0 0.0 UD M 2 2.0 0.0 ED M 3 2.0 1.0 ED

4 Results

In Fig. 5 , we present the model M 1 for a constant surface density disc (UD). The density contour and the velocity vectors for the wind are shown in four snapshots in Fig. 5 up to a time t = 98 Myr (corresponding to t = 2 in computational time units). There are a few aspects of the gaseous flow that we should note here. First, the disc and the outflowing gas in this case have no rotation ( vφ = 0). In the absence of the centrifugal force due to rotation, which might have reduced the radial gravitational force, there is a net radial force driving the gas inwards. At the same time, the radiation force, here characterized by Γ 0 = 2, propels the gas upwards (the radial component of radiation being weak). The net result after a few Myr is that the gas in the region near the pole moves in the positive z -direction, and there is a density enhancement inside a cone around the pole, away from which the density and velocities decrease.

Figure 5

M 1 : logarithmic density contours for a radiation-driven wind from UD for four snapshots running up to t = 98 Myr, with velocity vectors shown with arrows. Densities are colour-coded according to the computational unit of density, .

Figure 5

M 1 : logarithmic density contours for a radiation-driven wind from UD for four snapshots running up to t = 98 Myr, with velocity vectors shown with arrows. Densities are colour-coded according to the computational unit of density, .

Also, because of the strong gravity of the bulge, the gas tends to get trapped inside the bulge region, and even the gas at larger r tends to get dragged towards the axis. This region puffs due to the accumulation of matter. Ultimately, the radiative force drives matter outwards in the form of a plume.

Next, we change the disc mass distribution and simulate the case of wind driven out of an ED. We show the results in Fig. 6 . Since both gravity and radiation forces in this case of ED are quasi-spherical in nature, in the final snapshot the flow appears to follow almost radial streamlines. Although in the vicinity of the disc, the injected matter still falls towards the axis, but this is not seen at large height as was seen in the previous case of M 1 . This makes the wind cone of rising gas more diverging than in the case of UD (M 1 ).

Figure 6

M 2 : logarithmic density contours for a radiation-driven wind from ED for four snapshots running up to t = 98 Myr, with velocity vectors shown with arrows.

Figure 6

M 2 : logarithmic density contours for a radiation-driven wind from ED for four snapshots running up to t = 98 Myr, with velocity vectors shown with arrows.

4.1 Rotating wind from exponential disc

The direction of the fluid flow in M 1 and M 2 is by and large towards the axis, and this flow is mitigated in the presence of rotation in the disc and injected gas. In the next model M 3 , we consider rotating matter being injected into the computational domain and which follows a vφ distribution given by equation (23) . This is reasonable to assume since the disc from which the wind is supposed to blow is itself rotating. In M 3 , we simulate rotating gas being injected above an ED and being driven by a radiation force of Γ 0 = 2. We present nine snapshots of the M 3 case in Fig. 7 .

Figure 7

M 3 : contours of log 10 (ρ) and v -field of a radiation-driven wind with Γ 0 = 2.0 from an ED. t = 2 corresponds to 98 Myr.

Figure 7

M 3 : contours of log 10 (ρ) and v -field of a radiation-driven wind with Γ 0 = 2.0 from an ED. t = 2 corresponds to 98 Myr.

The first six snapshots of Fig. 7 show the essential dynamics of the outflowing gas. The fast rotating matter from the outer disc is driven outwards because the radial gravity component is balanced by rotation. Near the central region, the rotation and also the radial force components are small. Therefore, the gas is mostly driven vertically. The injected gas reaches a vertical height of ∼5 kpc in a time-scale of ∼37 Myr ( t = 0.75). The flow reaches a steady state after ∼60 Myr ( t = 1.25). In the steady state we find a rotating and mildly divergent wind.

We show the azimuthal velocity contours in Fig. 8 in colour for the fully developed wind (last snapshot in M 3 ), and superpose on it the contour lines of ρ. The density contours clearly show a conical structure for outflowing gas. The rotation speed of the gas peaks at the periphery of the cone, and is of the order of ∼50–100 km s −1 . Compared to the disc rotation speed, the rotation speed of the wind region is somewhat smaller. In other words, we find the wind mostly consisting of low angular momentum gas lifted from the disc.

Figure 8

The rotation velocity vφ for the case M 3 at a time of 98 Myr is shown in colours. Contour lines of log 10 (ρ) are plotted over it.

Figure 8

The rotation velocity vφ for the case M 3 at a time of 98 Myr is shown in colours. Contour lines of log 10 (ρ) are plotted over it.

We plot the velocity of the gas close to the axis in Fig. 9 for different times in this model (M 3 ), using . The velocity profile in the snapshots at earlier time fluctuates at different height, but becomes steady after t ≥ 1.5, as does the density profile.

Figure 9

The axial velocity vz (0 + , z ) with z at different time-steps for the model M 3 . t = 2.0 corresponds to a time of 98 Myr.

Figure 9

The axial velocity vz (0 + , z ) with z at different time-steps for the model M 3 . t = 2.0 corresponds to a time of 98 Myr.

We have run this particular case of ED with rotation (model M 3 ) for different values of Γ 0 . In order to illustrate the results of these runs, we plot the z -component of velocity [ vz (0 + , 10 kpc)] at 10 kpc and at simulation time, t = 2 as a function of Γ 0 in Fig. 10 . We find that significant wind velocities are obtained for Γ 0 ≳ 1.5 and wind velocities appear to rise linearly with Γ 0 after this critical value is achieved. Sharma et al. (2011) found this critical value to be Γ 0 ∼ 2 for a constant density disc and wind launched above the bulge. For the realistic case of an ED, we find in the present simulation the critical value to be somewhat smaller than but close to the analytical result. The important point is that the critical Γ 0 is not unity. This is because the parameter Γ 0 is not a true Eddington parameter since it is defined in terms of disc gravity and radiation, whereas halo and bulge also contribute to gravity.

Figure 10

The axial velocity vz (0 + , 10 kpc) in simulation units vref = 200 km s −1 with Γ 0 , at a time t ∼ 10 2 Myr.

Figure 10

The axial velocity vz (0 + , 10 kpc) in simulation units vref = 200 km s −1 with Γ 0 , at a time t ∼ 10 2 Myr.

5 Discussions

Our simulation differs from earlier works (e.g. Suchkov et al. 1994 ) mainly in that we specifically target warm outflows and the driving force is radiation pressure. Most of the previous simulations of galactic wind have used energy injected from SNe blasts as a driving force. However, with the ideas presented in Murray et al. (2005) , which worked out the case of radiation pressure in a spherical symmetric set-up, it becomes important to study the physics of this model in an axisymmetric set-up, as has been done analytically by Sharma et al. (2011) (see also Zhang & Thompson 2010 ). We have also tried to capture all features of a typical disc galaxy like a bulge and a dark matter halo, and a rotating disc. Recent analytical works ( Sharma & Nath 2011 ) and simulations ( Hopkins et al. 2012 ) have shown that outflows from massive galaxies ( Mhalo ≥ 10 12 M ) have different characteristics than those from low-mass galaxies. Outflows from massive galaxies are mostly driven by radiation pressure, and the fraction of cold gas in the haloes of massive galaxies is large ( van de Voort & Schaye 2011 ). Our simulations presented here address these outflows in particular.

We have parametrized our simulation runs with the disc Eddington factor Γ 0 , and it is important to know the corresponding luminosity for a typical disc galaxy, or the equivalent SFR. For a typical opacity of a dust and gas mixture (κ∼ 200 cm 2 g −1 ) ( Draine 2011 ), the corresponding mass-to-light ratio requirement for Γ 0 ≳ 1.5 is that M / L ≤ 0.03. Sharma et al. (2011) showed that for the case of an instantaneous star formation, Γ 0 ≳ 2 is possible for an initial period of ∼10 Myr after the starburst. However, for a continuous star formation, which is more realistic for disc galaxies, Sharma & Nath (2012) found that only ULIGs, with SFR larger than ∼100 M yr −1 and which are also massive, are suitable candidates for such large values of Γ 0 , and for radiatively driven winds.

The results presented in the previous sections show that the outflowing gas within the central region of a few kpc tends to stay close to the pole, and does not move outwards because of its low angular momentum. This makes the outflow somewhat collimated. Although outflows driven by SN-heated hot wind also produce a conical structure (e.g. Fragile et al. 2004 ) emanating from a breakout point of the SN remnants, there is a qualitative difference between this case and that of radiatively driven winds as presented in our simulations. While it is the pressure of the hot gas that expands gradually as it comes out of a stratified atmosphere, in the case of a radiation-driven wind, it is the combination of mostly the lack of rotation and almost vertical radiation driving force in the central region that produces the collimation effect.

We also note that the conical structure of rotation in the outflowing gas is similar to the case of outflow in M82 ( Greve 2004 ), where one observes a diverging and rotating periphery of conical outflow.

We have not considered radiative cooling in our simulations, since for typical density in the wind the radiative cooling time is shorter or comparable than the dynamical time. However, there are regions of higher density close to the base, and radiative cooling can be important there. Also radiation transfer may be important, particularly in optically thick environments of ULIGs and radiation flux may get reduced with distance ( Novak, Ostriker & Ciotti 2012 ). We hope to address these issues in a future paper.

From our results of the exponential and rotating disc model, we find the wind comprising low angular momentum gas lifted from the disc. It is interesting to note that recent simulations of SNe-driven winds have also claimed a similar result ( Governato et al. 2010 ). Such loss of low angular momentum gas from the disc may have important implication for the formation and evolution of the bulge, since the bulge population is deficient in stars with low specific angular momentum. Binney, Gerhardt & Silk (2001) have speculated that outflows from disc that preferentially removes low angular momentum material may resolve some discrepancies between observed properties of disc and results of numerical simulations.

As a caveat, we should finally note that the scope and predictions of our simulation are limited by the simple model of disc radiation adopted here. In reality, radiation from discs is likely to be confined in the vicinity of star clusters, and not spread throughout the disc as we have assumed here. This is likely to increase the efficacy of radiation pressure, which is not possible within the scope of an axisymmetric simulation.

6 Summary

We have presented the results of hydrodynamical (Eulerian) simulations of radiation-driven winds from disc galaxies. After studying the cases of winds from a constant surface density disc and ED without rotation, we have studied a rotating outflow originating from an ED with rotation. We find that the outflow speed increases rapidly with the disc Eddington parameter Γ 0I /(2 cG Σ) for Γ 0 ≥ 1.5, consistent with theoretical expectations. The density structure of the outflow has a conical appearance, and most of the outflowing gas consists of low angular momentum gas.

We thank Yuri Shchekinov for constructive comments and critical reading of the manuscript. IC acknowledges the hospitality of the Astronomy and Astrophysics Group of Raman Research Institute, where the present work was conceived. DR was supported by National Research Foundation of Korea through grant 2007-0093860.

REFERENCES

Bianchi
S.
Ferrara
A.
,
2005
,
MNRAS
,
358
,
379
Binney
J.
Gerhardt
O.
Silk
J.
,
2001
,
MNRAS
,
321
,
471
Burke
A. J.
,
1968
,
MNRAS
,
140
,
241
I.
,
2005
,
MNRAS
,
356
,
145
Chevalier
R. A.
Clegg
A. W.
,
1985
,
Nat
,
317
,
44
Chiao
R. Y.
Wickramasinghe
N. C.
,
1972
,
MNRAS
,
159
,
361
Cooper
J. L.
Bicknell
G. V.
Sutherland
R. S.
Bland-Hawthorn
J.
,
2008
,
ApJ
,
674
,
157
Davies
J. I.
Alton
P.
Bianchi
S.
Trewhella
M.
,
1998
,
MNRAS
,
300
,
1006
Dekel
A.
Silk
J.
,
1986
,
ApJ
,
303
,
39
Dickey
J. M.
Lockman
F. J.
,
1990
,
ARA&A
,
28
,
215
Draine
B. T.
,
2011
,
ApJ
,
732
,
100
Dubois
Y.
Teyssier
R.
,
2008
,
A&A
,
477
,
79
Fragile
P. C.
Murray
S. D.
Lin
D. N. C.
,
2004
,
ApJ
,
617
,
1077
Fujita
A.
Mac Low
M.
Ferrara
A.
Meiksin
A.
,
2004
,
ApJ
,
613
,
159
Fujita
A.
Martin
C. L.
Mac Low
M.
Kimberly
C. B.
Weaver
R.
,
2009
,
ApJ
,
698
,
693
Governato
F.
et al. ,
2010
,
Nat
,
463
,
203
Greve
A.
,
2004
,
A&A
,
416
,
67
Harten
A.
,
1983
,
J. Comput. Phys.
,
49
,
357
Heckman
T. M.
,
2002
, in
Mulchaey
J. S.
Stocke
J. T.
, eds,
ASP Conf. Ser. Vol. 254. Extragalactic Gas at Low Redshift
.
Astron. Soc. Pac.
, San Francisco, p.
292
Hopkins
P. F.
Quataert
E.
Murray
N.
,
May 2012
, MNRAS, 421, 3522
Johnson
H. E.
Axford
J. I.
,
1971
,
ApJ
,
165
,
381
Kang
H.
Ostriker
J. P.
Cen
R.
Ryu
D.
Hernquist
L.
Evrard
A. E.
Bryan
G. L.
Norman
M. L.
,
1994
,
ApJ
,
430
,
83
Kohji
T.
Ikeuchi
S.
,
1988
,
ApJ
,
330
,
695
Larson
R. B.
,
1974
,
MNRAS
,
169
,
229
Lynds
C.
Sandage
A.
,
1963
,
ApJ
,
137
,
1005
Mac Low
M.
Ferrara
A.
,
1999
,
ApJ
,
513
,
142
Martin
C. L.
,
2005
,
ApJ
,
621
,
227
Mathews
W. G.
Baker
R. G.
,
1971
,
ApJ
,
170
,
241
Mo
H. J.
Mao
S.
White
S. D. M.
,
1998
,
MNRAS
,
295
,
319
(MMW98)
Molteni
D.
Ryu
D.
Chakrabarti
S. K.
,
1996
,
ApJ
,
470
,
460
Murray
N.
Quataert
Q.
Thompson
T. A.
,
2005
,
ApJ
,
618
,
569
Murray
N.
Ménard
B.
Thompson
T. A.
,
2011
,
ApJ
,
735
,
66
Nath
B. B.
Silk
J.
,
2009
,
MNRAS
,
396
,
L90
Navarro
J. F.
Frenk
C. S.
White
S. D. M.
,
1997
,
ApJ
,
490
,
493
Novak
G. S.
Ostriker
J. P.
Ciotti
L.
,
May 2012
,
preprint (arXiv:1203.6062)

Oppenheimer
B. D.
Davé
R.
,
2006
,
MNRAS
,
373
,
1265
Rupke
D. S.
Veilleux
S.
Sanders
D. B.
,
2002
,
ApJ
,
570
,
588
Rupke
D. S.
Veilleux
S.
Sanders
D. B.
,
2005
,
ApJS
,
160
,
115
Ryu
D.
Ostriker
J. P.
Kang
H.
Cen
R.
,
1993
,
ApJ
,
414
,
1
Ryu
D.
Brown
G. L.
Ostriker
J. P.
Loeb
A.
,
1995
,
ApJ
,
452
,
364
Savage
B. D.
Sembach
K. R.
Lu
L.
,
1997
,
AJ
,
113
,
2158
Seaquist
E. R.
Clark
J.
,
2001
,
ApJ
,
552
,
133
Shapley
A. E.
Steidel
C. C.
Pettini
M.
K. L.
,
2003
,
ApJ
,
588
,
65
Sharma
M.
Nath
B. B.
,
May 2012
, ApJ, 750, 55
Sharma
M.
Nath
B. B.
Shchekinov
Y.
,
2011
,
ApJ
,
736
,
L27
Sofue
Y.
Reuter
H. P.
Krause
M.
Wielebinski
R.
Nakai
N.
,
1992
,
ApJ
,
395
,
126
Strickland
D. K.
Stevens
I. R.
,
2000
,
ApJ
,
314
,
511
Strickland
D. K.
Heckman
T. M.
Colbert
E. J. M.
Hoopes
C. G.
Weaver
K. A.
,
2004
,
ApJS
,
151
,
153
Suchkov
A. A.
Berman
V. G.
Heckman
T. M.
Balsara
D. S.
,
1994
,
ApJ
,
430
,
511
Suchkov
A. A.
Berman
V. G.
Heckman
T. M.
Balsara
D. S.
,
1996
,
ApJ
,
463
,
528
Sutherland
R. M.
Dopita
M. A.
,
1993
,
ApJS
,
88
,
253
Tielens
A. G. G. M.
McKee
C. F.
Seab
C. G.
Hollenbach
D. J.
,
1994
,
ApJ
,
431
,
321
Tomisaka
K.
Bregman
J. N.
,
1993
,
PASJ
,
45
,
513
van de Voort
F.
Schaye
J.
,
2011
, preprint (arXiv:1111.5039)
Veilleux
S.
Cecil
G.
Bland-Hawthorn
J.
,
2005
,
ARA&A
,
43
,
769
Veilleux
S.
Rupke
D. S. N.
Swaters
R.
,
2009
,
ApJ
,
700
,
L149
Walter
F.
Weiss
A.
Scoville
N.
,
2002
,
ApJ
,
580
,
L21
Wesmoquette
et al. ,
2009
,
ApJ
,
696
,
192
West moquette
M. S.
Smith
L. J.
Gallacher
J. S.
Trancho
G.
Bastian
N.
Konstantopoulos
I. S.
,
2009
, ApJ, 696, 192
Zhang
D.
Thompson
T. A.
,
2010
, preprint(arXiv:1005.4691)

Appendix

APPENDIX A:  COMPUTATION OF DISC GRAVITY

Consider a razor-thin disc in r φ plane as illustrated in Fig. A1 . Now our task is to calculate the force components at any arbitrary point above the disc. Let us consider an annulus of the disc between r ′ and . The area of the element at point P ( r ′, φ′, 0) is . Also take a field point Q ( r , 0, z ) above the disc plane. The azimuthal coordinate of Q is taken to be zero for simplicity as we know that azimuthal force components are zero due to symmetry. Let QN and QM be perpendiculars from Q on the x - and z -axis, respectively. So we can write

(A1)
The gravitational force due to the small area element at P is given by
(A2)

Figure A1

Schematic diagram for the calculation of gravitational force due to disc in the xy -plane. We consider an annulus in the disc and an element of area around the point P ( r ′, φ, 0) in this annulus is considered here in order to compute the force at a point Q ( r , 0, z ) whose azimuthal coordinate φ= 0. The point S ( r ′cos φ, 0, 0) is the foot of the perpendicular drawn from P on the x -axis. The point S′ ( r ′, 0, z ) is at the intersection of the vertical from S (along z -axis) and the line parallel to x -axis at height z . The angle , and .

Figure A1

Schematic diagram for the calculation of gravitational force due to disc in the xy -plane. We consider an annulus in the disc and an element of area around the point P ( r ′, φ, 0) in this annulus is considered here in order to compute the force at a point Q ( r , 0, z ) whose azimuthal coordinate φ= 0. The point S ( r ′cos φ, 0, 0) is the foot of the perpendicular drawn from P on the x -axis. The point S′ ( r ′, 0, z ) is at the intersection of the vertical from S (along z -axis) and the line parallel to x -axis at height z . The angle , and .

Here Σ( r ′) is the surface density of the disc. Now the z -component of this force is

(A3)
To calculate the radial component, let PS be the perpendicular from P on the x -axis. Then, we have sin∠SPN = SN/PN = ( )/PN. Component of the force along the direction of PN is
(A4)