© ESO, 2014

1. Introduction

Originating in the Galactic centre and extending up to a distance of ≈15 kpc from the Galactic plane, there is a large bipolar structure with counterparts in different wavelengths, from radio to gamma rays. This structure was seen for the first time by ROSAT in X-rays (Snowden et al. 1997). An overlapping extended source was also detected in radio by WMAP (Finkbeiner 2004), and a few years later Fermidetected a somewhat smaller structure in the GeV range, the Fermi bubbles, with a height of ≈8 kpc (Su et al. 2010). A hint of an extension of the GeV emission beyond the Fermi bubbles was also found on scales similar to those of the ROSAT emission (Ackermann et al. 2012). S-PASS (Carretti et al. 2013) and WMAP (Jones et al. 2012) detected a large-scale, partly ordered magnetic-field structure placed at 20°−50° from the Galactic plane in the Fermibubble region, which is larger than the bubbles themselves.

The origin of the Fermi bubbles is still under debate. There are two main hypotheses, which were already discussed in Su et al. (2010): active galactic nucleus (AGN) activity, or a bipolar galactic wind. In fact, a combined origin of central accretion activity and a galactic wind fed by star formation may be possible, as discussed in Carretti et al. (2013). In that work, the authors adopted a leptonic scenario for the Fermi bubbles emission, whereas Crocker & Aharonian (2011) considered hadronic radiation.

Numerical simulations of the dynamics of the Fermi bubbles can be useful to unveil the relevant physical processes underlying the detected radiation. In this regard, numerical hydrodynamical simulations of the Fermi bubbles were performed in 2D by Guo & Mathews (2012) and Guo et al. (2012) and 3D by Yang et al. (2012). In the latter paper, the authors took into account anisotropic diffusion of cosmic rays to explain the sharp edges of the Fermi bubbles in gamma rays. In Guo & Mathews (2012), the X-rays detected by ROSAT were explained as thermal radiation from the surrounding medium heated by a shock, and the gamma rays from the Fermi bubbles as coming from relativistic electrons in the ejecta.

Motivated in particular by the recent detection of a partly-ordered magnetic-field structure in the Fermibubble region, we here present the results of 3D magnetohydrodynamical (MHD) simulations of the structures associated with the Fermi bubbles in the framework of a mini-AGN event in which material is ejected from the central Galactic black hole (see Bland-Hawthorn et al. 2013, for recent observational support for the AGN scenario). We interpret in this context the recent observational findings in radio, in particular concerning the large-scale magnetic field, and suggest a non-thermal nature for the X-ray radiation detected by ROSAT, which is produced by electrons accelerated in the shock driven by the ejecta in the surrounding medium. The GeV emission from the Fermi bubbles might result from particle acceleration in the shock wave of a second explosive event in a recurrent-activity scenario, or particle acceleration within the ejecta, in which shocks, turbulence, and suitable conditions for magnetic reconnection are present. Throughout, we adopted the notation A_b = A/10^b, where A has cgs units.

2. Physical scenario

We present an MHD model for the formation of the large-scale magnetic-field structure in the Fermibubble region in the context of accretion-driven explosive events. The evolution and observational properties of the Fermi bubbles and accompanying structures can be strongly influenced by the presence of this large-scale magnetic field.

We assumed that the very central region of the Galaxy is surrounded by a rotating gaseous torus or disc with outer radius R_D, effective thickness or inner radius of the torus r_D, and rotation speed υ_D. Loops of magnetic field thread the central region. These magnetic loops, of strength B_l, are anchored in the disc and surround the central region, as shown in Fig. 1.

Fig. 1 Sketch of the initial configuration: the opaque horizontal torus is a gaseous disc with the distance from the centre of the tube to the centre of the torus RD = 1 kpc and the radius of the tube rD = 0.2 kpc; the translucent vertical ring is a magnetic loop with the same radius and thickness as that of the gaseous disc. The central sphere of radius rs = 0.2 kpc represents a central region with high pressure and velocity: the ejection region.

In this scenario, the formation of the Fermibubble region proceeds as follows: a molecular disc, with a lifetime ~10⁷ years (Molinari et al. 2011), surrounds the Galactic centre. A magnetic-field loop anchored in the disc also crosses the disc polar regions. At a given time, and possibly episodically, the accretion rate in the central black hole increases dramatically because of for instance the capture of a molecular cloud. Such an accretion event can launch a powerful outflow in the polar direction, as occurs in AGN. The ejecta will stretch the magnetic-field loops that cross the polar regions, and push the external medium up- and sidewards, shocking it. Simultaneously, the external layers of the whole structure, where the ejecta and the external medium are in contact, will be twisted by the rotating gaseous disc through the surrounding magnetic field anchored in the disc and will also be sheared by differential rotation. As shown below, this scenario can produce the partly ordered magnetic-field structure similar to that found in S-PASS.

Fig. 2 Distribution of the gas velocity modulus and the magnetic-field lines for two models: on the left panel, Etot = 3 × 1055 erg at a time t = 2.3 Myr; on the right panel, Etot = 1055 erg at t = 3.4 Myr. The velocity is presented in units of 108 cm/s and the magnetic-field strength (shown in colour) in units of 400 μG.

3. Simulations

The simulations presented here were implemented in 3D Cartesian geometry using the PLUTO code (Mignone et al. 2007, 2012), the piece-parabolic method (PPM, Colella & Woodward 1984; Mignone et al. 2005), an HLLD Riemann solver (Miyoshi & Kusano 2008), and applying AMR using the Chombo code (Colella et al. 2009). The flow was approximated as an ideal, adiabatic gas with a magnetic field, one-particle species, and a polytropic index of 5/3. The adopted resolution was 48 × 48 × 96 cells and three levels of AMR, which gives an effective resolution of 384 × 384 × 768 cells. The size of the domain was x ∈ [0,18 kpc], y ∈ [0,18 kpc], and z ∈ [0,36 kpc]. The calculations were carried out in the Moscow State University cluster Chebyshev. The visualization of the results were arranged with VisIt.

3.1. Initial setup

In the centre of our galaxy there is a massive gaseous disc of ~3 × 10⁷ M_⊙, a radius of 100 pc, and a rotation speed ~80 km s^-1 (see Molinari et al. 2011). The initial configuration of the simulation sets an horizontal gaseous torus in the plane XY (the disc), and a vertical magnetic-field loop in the plane XZ. The torus radius and effective thickness were taken to be R_D = 1 kpc and r_D = 0.2 kpc, respectively. We fix the central particle density of the torus to 10³ cm^-3, and assumed constant angular momentum and hydrostatic equilibrium to determine the torus pressure and temperature distribution. We used a truncated Keplerian gravitational potential φ ∝ 1/max(r,0.6 kpc). This potential will not significantly influence the calculation results because the speed of the shock wave at all radii is significantly higher than the escape velocity. The magnetic-loop radius and thickness were taken to be equal to those of the torus (see Fig. 1). The initial strength of the magnetic field in the loop was 100 μG, similar to the observed value in the Galactic centre (Crocker et al. 2010, 2011). Numerical limitations led us to adopt a larger size of the disc/torus and the associated magnetic loop. This also implied that the disc was simulated to be rotating in the gravitational potential of the central region with a rather high speed, υ_D = 800 km s^-1, to obtain a disc angular velocity equal to the actual angular velocity of the molecular disc in the centre of the Galaxy. Hence, in practice, the main difference between our simulation and a more realistic setup is the width of the inner ejection region, which is expected to have a minor impact on the long-term evolution of the simulated structures. In addition, the somewhat unrealistic torus properties do not allow a precise, quantitative comparison between the simulated magnetic-field strength and the observed values. However, given that the magnetic field is not dynamically relevant, it is determined by the fluid evolution, and the obtained global field geometry is expected to be reliable. The lifetime of the central engine activity was assumed to be much shorter than that of the disc or the dynamical scale of the simulation.

The explosive activity produces an axial outflow from a central sphere of radius r_s = 0.2 kpc. The outflow, initially confined to this sphere, was set with two different total energies for one simulation with one ejection: E_tot = 10⁵⁵ and 3 × 10⁵⁵ erg; a particle density of 1 cm^-3, and a velocity of ${\mathit{v}}_{\mathrm{e}}\mathrm{=}\mathrm{3.5}\hspace{0.17em}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}\mathrm{sign}\mathrm{\left(}\mathrm{z}\mathrm{\right)}{\mathit{E}}_{\mathrm{56}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ cm s^-1 at ejection, which is directed along the Z-axis. The gas density of the surrounding medium was set as in Miller & Bregman (2013): ${\mathit{n}}_{\mathrm{bg}}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{0}}}{{\left[\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{(}\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{C}}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{(}\mathit{z}\mathit{/}{\mathit{z}}_{\mathrm{C}}{\mathrm{)}}^{\mathrm{2}}\right]}^{\mathrm{3}\mathit{\beta }\mathit{/}\mathrm{2}}}\mathit{,}$ where n₀ = 0.46 cm^-3, R is the cylindrical radius, z the distance from the Galactic plane, and β = 0.71, R_C = 0.42 kpc and z_C = 0.26 pc are normalizing constants. The temperature of this medium was fixed to 1.2 × 10⁶ K. Assuming that accretion events may be recurrent, we performed an additional simulation with two active episodes. An amount of energy E_tot = 10⁵⁵ erg was now injected twice, at t = 0 and at t = 2.5 Myr. The central engine activity was taken to be shorter than the flow-crossing time, r_s/v_e ~ 10⁴ yr. Because they are much shorter than any relevant dynamical process in the simulation, the activity periods can indeed be considered as discrete events. Note that the formation timescale of this magnetic field structure is ~100 pc / 80 km s^-1 ≈ 1 Myr, which implies that the magnetic loop can continue to be generated between events.

3.2. Results

The distributions of velocity and magnetic-field lines are presented in Fig. 2 for two models, E_tot = 3 × 10⁵⁵ and E_tot = 10⁵⁵ erg, at the time when the size along the Z-direction of the whole structure reached ≈15 kpc. The expansion of the high-pressure region inflates the bubbles and stretches the magnetic-field lines, and meanwhile the rotation of the magnetic-line footpoints twist the magnetic field into a spiral structure. This spiral structure is similar to that observed in Carretti et al. (2013).

The ejected matter from the central engine forms a bipolar, elongated and wide structure. The shock wave travelling through the Galactic halo forms an external bubble of shocked external medium almost as wide as it is long. Surrounded by this shocked external medium, the ejected matter forms a narrower second bubble with its neck at the Galactic centre (see Fig. 3). Covering the inner bubble, at the contact layer between the ejecta and the shocked external medium, the pushed magnetic field is threaded, sheared and anchored in the rotating disc. The magnetic structure induces minor perturbations at the contact layer. In the simulation with two active episodes a second bubble forms inside the first one, as shown in Fig. 4. In this case, the sheared magnetic lines anchored in the molecular disc surround the shock wave formed by the second bubble on the earlier one. An interesting property of the computed solution is the higher speed at the top of the bubble shocks than at their sides. This effect is produced by a decreasing background-density profile.

4. Discussion

Fig. 3 Distribution of the ejected matter tracer (opaque blue-green) and pressure surface (transparent magenta) indicating the position of the shock wave from the model with Etot = 3 × 1055 erg at t = 2.3 Myr.

Fig. 4 Distribution of the gas pressure (bottom and top) and magnetic-field lines (top) at t = 3.4 Myr for a model with two active episodes, with Etot = 1055 erg each at times t = 0 and t = 2.5 Myr. The pressure is presented in units of 1.67 × 10-8 pa and the magnetic field (shown in colour) in units of 400 μG.

We proposed an accretion-ejection, explosive scenario for the formation of the large-scale magnetic-field structure around the Fermi bubbles. As shown in Fig. 3, one explosion produces very extended ejecta, much larger than the Fermi bubbles themselves, whose head reaches up to the external boundary of the whole structure. This mismatch in size may be explained by invoking acceleration and confinement of the GeV-emitting particles only deep into the ejecta. Another possibility would come from accretion-ejection events recurring on timescales ≳10⁶ Myr. Such a timescale, and the energetics derived here, would imply an average injection luminosity of ~10⁴¹ erg s^-1, compatible with cosmic-ray injection. As shown in Fig. 4, the two-episodes scenario, with this characteristic timescale and energy budget, nicely reproduces the overall morphology of the whole Fermibubble region, in particular the sharp edge of the Fermi bubbles themselves¹. It is worth noting that independently of the simulation, the assumption of a magnetic structure anchored in the central molecular disc already provides a constraint on the age of the Fermi bubbles of ≲10⁷ yr. If this were true, such a short lifetime would disfavour the hadronic scenario. This is so because of the long cooling time of proton-proton collisions (Kelner et al. 2006): ${\mathit{t}}_{\mathrm{pp}}\mathrm{=}{\mathrm{10}}^{\mathrm{15}}\mathit{/}\mathit{n}\mathit{>}{\mathrm{10}}^{\mathrm{18}}{\mathit{n}}_{-4}^{-1}\mathrm{s}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{10}}{\mathit{n}}_{-4}^{-1}$ yr, ~10⁴ times longer than the dynamical time of the magnetic field. This low efficiency would require a total energy only in high-energy protons ~L_γt_pp/0.17 = 2 × 10⁵⁶ ergs, ten times higher than the total energetics of the Fermi bubbles if they are related to the large-scale magnetic-field structure, as assumed here. Accounting for dynamical constraints, that is, to explain the growth of the structure in the surrounding medium on the timescales considered here, one derives an estimate on the luminosity injected into the Fermi bubbles of ~10⁴¹ ergs. This is ~10³ times higher than the energy required to explain the gamma-ray luminosity, which means that the energy budget is not tight when considering a leptonic origin of the gamma-ray emission.

We now briefly consider a framework for the broadband detected radiation. Electrons and protons can be accelerated at the shock produced by the ejected matter in the surrounding medium. As just argued, protons are disfavoured by the relatively short timescales involved in the adopted scenario. For electrons, on the other hand, this is not a problem, with the most efficient radiation process at low energies being synchrotron, and at high energies, inverse Compton (IC) with soft Galactic photons of ~1 eV. The maximum energy of the photons produced via synchrotron from shock-accelerated electrons can be estimated as ϵ_s ≈ 0.12(υ_s,8)² η^-1 keV (Aharonian & Atoyan 1999), where υ_s is the forward-shock speed. However, given the slow steepening with energy above ϵ_s of the synchrotron spectrum of shock-accelerated electrons (Zirakashvili & Aharonian 2007), the effective maximum energy of the synchrotron photons can be estimated to be about 10 × ϵ_s. With a speed of the shock wave of ~2 × 10⁸ cm s^-1, synchrotron X-rays can reach up to few keV. This synchrotron radiation might be able to explain the X-ray shell or arch found by ROSAT (Snowden et al. 1997), which would then have a non-thermal origin instead of a thermal one (both possibilities were discussed for instance); IC emission from the same electrons might explain the hint of GeV radiation found on similar scales (see Fig. 6 in Ackermann et al. 2012). A diffuse lower-energy electron population, extended also beyond the Fermi bubbles, but partially embedded in the spiral-like magnetic lines, might explain the WMAP extended source. The Fermi bubbles themselves would have an IC origin. They would remain confined to the fresher material of the most recent episode, or alternatively, to deeper regions of the whole region if only one explosion took place, as was discussed by Guo & Mathews (2012). A hardening of the GeV radiation along the Z-direction has been observed for the Fermi bubbles (Hooper & Slatyer 2013). This hardening might be explained by the varying speed of the bubble shocks mentioned in Sect. 3.2, affecting (at least slightly) the slope of the distribution of the GeV emitting particles. This hardening might also be consistent with a stronger shock at the top of the bubbles than at their sides.

Acknowledgments

We thank the anonymous referee for useful and constructive comments. The calculations were carried out in the cluster of Moscow State University Chebyshev. We thank Andrea Mignone and the PLUTO team for the possibility to use the PLUTO code. We also thank the Chombo team for the possibility to use the Chombo code. We would like to thank the VisIt visualization packet team. B.M.V. acknowledges partial support by RFBR grant 12-02-01336-a. V.B-R. acknowledges support by DGI of the Spanish Ministerio de Economía y Competitividad (MINECO) under grants AYA2010-21782-C03-01 and FPA2010-22056-C06-02. V.B-R. acknowledges financial support from MINECO through a Ramón y Cajal fellowship. This research has been supported by the Marie Curie Career Integration Grant 321520.

References

All Figures

Fig. 1 Sketch of the initial configuration: the opaque horizontal torus is a gaseous disc with the distance from the centre of the tube to the centre of the torus RD = 1 kpc and the radius of the tube rD = 0.2 kpc; the translucent vertical ring is a magnetic loop with the same radius and thickness as that of the gaseous disc. The central sphere of radius rs = 0.2 kpc represents a central region with high pressure and velocity: the ejection region.

In the text

	Fig. 2 Distribution of the gas velocity modulus and the magnetic-field lines for two models: on the left panel, Etot = 3 × 1055 erg at a time t = 2.3 Myr; on the right panel, Etot = 1055 erg at t = 3.4 Myr. The velocity is presented in units of 108 cm/s and the magnetic-field strength (shown in colour) in units of 400 μG.
In the text

	Fig. 3 Distribution of the ejected matter tracer (opaque blue-green) and pressure surface (transparent magenta) indicating the position of the shock wave from the model with Etot = 3 × 1055 erg at t = 2.3 Myr.
In the text

	Fig. 4 Distribution of the gas pressure (bottom and top) and magnetic-field lines (top) at t = 3.4 Myr for a model with two active episodes, with Etot = 1055 erg each at times t = 0 and t = 2.5 Myr. The pressure is presented in units of 1.67 × 10-8 pa and the magnetic field (shown in colour) in units of 400 μG.
In the text