The morphology of the middle-aged supernova remnant, Cygnus Loop, seen in X-rays, is peculiar, with a blowout in the south region and other irregular features, such as a bump in the west, a limb with a planar morphology in the east and asymmetry between the east and the west shock profiles of the blowout. The detailed process of the formation of the peculiar profile of the shock is still unclear. We perform three-dimensional hydrodynamical simulations for the remnant to revisit its evolution. In the simulations, the progenitor ejects an anisotropic, latitude-dependent wind, and travels in a direction that is not aligned with the symmetry axis of the wind. As a result, a cavity with a fringed structure is produced. The remnant has evolved in the cavity for about 104 yr. In the north-east, the shock has first encountered the bow shock, and this part corresponds to the bright north-eastern region. The south blowout is formed due to the shock travelling into the undisturbed wind, and the interaction of the shock with the cavity leads to the other peculiar features of the shock structure. The resulting profile of the remnant is consistent with that indicated in X-rays. It can be concluded that the supernova explosion occurred in the cavity produced by an anisotropic stellar wind experiencing two main phases with different wind velocities.

1 INTRODUCTION

The Cygnus supernova remnant (SNR) has an age of ∼104 yr and a large apparent size with a radius of ∼15 pc (Levenson et al. 1997, 1998) for a distance of 540 pc, which is derived from the detected proper motion of the Balmer-dominated filament in the north-east (Blair, Sankrit & Raymond 2005). Its morphology in X-rays is complicated, with a prominent limb in the north-east and some hollows along the shell (Levenson et al. 1998). Especially, a bright knot, called the XA region, exists on the east of the remnant (Hester & Cox 1986; McEntaffer & Brantseg 2011), and a blowout exists in the southern part of the rim (Levenson et al. 1997, 1998; Uchida et al. 2008). The origins of these irregularities are not well understood. The blowout region has been observed with XMM–Newton (Uchida et al. 2008), and the results show that the X-rays have two components. The low-kTe component comes from the interstellar medium, and the high-kTe component originates from the ejecta. Aschenbach & Leahy (1999) have suggested that the blowout region was produced as the shock wave travels into a lower-density medium, and the X-ray observations support this proposition.

According to the spectral analysis of the Cygnus Loop from the observations with XMM–Newton, the mass of the progenitor star is estimated to be ∼15 M (Tsunemi et al. 2007), and the remnant has evolved in the inhomogeneous medium (Zhou et al. 2010). The density distribution of the X-ray emitting plasma leads to the argument that the supernova has occurred within a bubble generated by the stellar wind of the progenitor star, and the shock has already encountered the wall of the cavity. However, the detailed formation of the morphology for the Cygnus Loop, with the blowout and the other irregularities evolved in such a wind-blown bubble, is still unresolved.

Numerical simulations have widely been adopted to study the dynamical evolution and the morphology of SNRs (Orlando et al. 2007, 2011), such as SN 1006 (Bocchino et al. 2011; Schneiter et al. 2015; Yu et al. 2015), RX J0852.0–4622 (Fang, Yu & Zhang 2014), G1.9+0.3 (Tsebrenko & Soker 2015) and Cassiopeia A (Orlando et al. 2016). For instance, Toledo-Roy et al. (2014a) investigated the dynamical properties and the X-ray shape for Kepler's SNR by assuming the remnant evolved in the wind bubble produced by the asymptotic giant branch (AGB) companion star. In their model, the wind of the companion star with a large runaway velocity is anisotropic, and the synthetic X-ray map can explain the observed X-ray morphology with a bright north-west shell and a bar-like structure. Moreover, Toledo-Roy et al. (2014b) numerically reproduced the peculiar radio shape of the SNR G352.7–0.1, which has a shell-like inner ring and an outer arc, by assuming that the progenitor exploded inside and near the border of a cloud. Recently, Meyer et al. (2015) investigated the asymmetry of SNRs evolving in an environment influenced by the wind of the progenitor star with a runaway velocity using two-dimensional hydrodynamical simulations. In the simulations, one side of the shock wave of the evolved remnant has collided with the bow shock, and that on the opposite direction expands into the undisturbed wind material.

In this paper, based on three-dimensional (3D) hydrodynamical simulations, we intend to explain the peculiar profile of the Cygnus Loop SNR. We assume that the progenitor wind, which is anisotropic with a latitude dependence, experiences two stages with different wind velocities during the evolution, and that the progenitor star has a spatial velocity misaligned with the polar axis of the wind. The organization of this paper is as follows. In Section 2, we describe the 3D model and the numerical simulations. In Section 3, we present the results of our simulations. Finally, in Section 4, we give a discussion and our conclusion.

2 NUMERICAL SET-UP

In this section, we present the numerical model for the evolutions of the anisotropic wind of the progenitor and the SNR within the wind-blown cavity. The 3D hydrodynamical numerical simulations are performed in a uniform grid based on the pluto code (Mignone et al. 2007, 2012).

2.1 Modelling the wind-blown cavity

To reproduce the asymmetric morphology of the remnant, the wind of the progenitor star with a velocity vs along the z-direction is assumed to be anisotropic with a latitude dependence, as in Raga et al. (2008), in which the assumption adopted is to model the AGB star Mira, and the polar axis of the star is inclined from the z-direction with an angle of α. The density decreases from the equator towards the pole with an equator-to-pole ratio ξ (Raga et al. 2008)
(1)
with a terminal velocity of
(2)
Here, vp is the terminal velocity of the wind at the pole, θ is the polar angle, and the function f is given by
(3)
The scaling constant A is
(4)
where | $\dot{M}_{\rm w}$ | is the mass-loss rate of the material blown from the star. The wind material expands freely to the interstellar medium with a uniform density of ρISM.

Two phases of the stellar wind of the progenitor are adopted to reproduce the peculiar shape of the remnant. The first phase lasts for a time of Tw1 with a terminal velocity at the pole of vp1. Then, the wind enters into the next stage, with a duration of Tw2 and a velocity of vp2, which is significantly larger than vp1.

2.2 Set-up of the ejecta

The dynamical evolution of the remnant is initiated by setting an ejecta that contains a kinetic energy Eej and a mass Mej in the simulation, with an inner core of a uniform density and an outer layer with a power-law distribution (Fraschetti et al. 2010), i.e.
(5)
The index s = 9 is used in this paper for the core-collapse supernova explosion. Moreover, the velocity of the matter in the ejecta at | $\boldsymbol {r}$ | is
(6)
and rc can be obtained by
(7)
Here, η = 3/7 is the mass ratio of the outer part of the ejecta to the inner part, and Rej = 1.0 pc is the initial radius of the ejecta. The density of the inner core is (Jun & Norman 1996; Fang & Zhang 2012)
(8)
and the velocity at the outer border of the ejecta v0 is given by
(9)

2.3 Modelling the SNR evolving in the cavity

The simulations are performed in a 3D Cartesian coordinate system. First, the wind of the progenitor star is simulated as material continually injected from a boundary with r = 1.0 pc in the two phases with a constant | $\dot{M}_{\rm w}$ | and durations of Tw1 and Tw2, respectively. Secondly, the density and the velocity distributions of the wind material are adopted as the initial conditions for the evolution of the ejecta.

The dynamics of the gas are solved with the Euler equations:
(10)
(11)
(12)
Here, P is the gas pressure and E is the total energy density,
t is the time, ρ = μmHn is the mass density, μ = 0.6 for the ionized gas with an assumption of a 10 : 1 H : He ratio, mH is the mass of a hydrogen atom, n is the gas number density, γ = 5/3 is the adiabatic index for non-relativistic gas, | $\boldsymbol {v}$ | is the gas velocity and Λ(T) is a tabulated cooling function (Dalgarno & McCray 1972; Toledo-Roy et al. 2009).

3 RESULTS

By assuming that the remnant evolves in the cavity produced by an anisotropic stellar wind from the progenitor star, hydrodynamical simulations are performed to reproduce the peculiar morphology of the Cygnus Loop observed in the X-ray band. Fig. 1 shows the X-ray morphology in the band 0.5–2.0 keV and the distribution of the hard/soft (0.5–2.0/0.1–0.3 keV) band ratio. Obviously, the X-ray images detected with ROSAT indicate a prominent feature that the interior of the Cygnus Loop is harder, and therefore hotter, than the edge. Besides the blowout in the south region, the limb of the Cygnus Loop in the X-ray images is approximately circular with some depressions on the edge. Especially, there is a prominent bump in the west, and the outer shock from the east to south-east is planar with a depression in the middle. Moreover, the shock to the east of the blowout is more circular with a much smaller curvature radius than that to the west.

The hard band (0.5–2.0 keV) mosaic of the Cygnus Loop (left panel) and the hard/soft band ratio (0.5–2.0/0.1–0.3 keV) of the remnant (right panel) obtained with ROSAT (Levenson, Graham & Snowden 1999). This image has been downloaded from http://imagine.gsfc.nasa.gov/Images/rosat/snr_cygloop.html. (Credit: S. L. Snowden)
Figure 1.

The hard band (0.5–2.0 keV) mosaic of the Cygnus Loop (left panel) and the hard/soft band ratio (0.5–2.0/0.1–0.3 keV) of the remnant (right panel) obtained with ROSAT (Levenson, Graham & Snowden 1999). This image has been downloaded from http://imagine.gsfc.nasa.gov/Images/rosat/snr_cygloop.html. (Credit: S. L. Snowden)

The simulations are initiated with a cubic Cartesian domain, which extends from −16 pc to 16 pc in each of the x-, y- and z-directions, and the environment flows in the −z-direction. The polar axis of the latitude-dependent wind lies on the yz-plane with an inclination angle of α to the z-axis. The evolutions of both the stellar wind and the remnant after the supernova explosion are calculated with 2563 grid cells. The circumstellar medium is homogeneous with a space velocity of vs towards the −z-direction, a density of nH and a temperature of T = 104 K. The stellar wind is injected in a sphere with a radius of 1 pc during a time Tw = 4.4 × 105 yr with a mass-loss rate of | $\dot{M}_{\rm w} = 10^{-5}$ | M yr−1. The matter distribution of the stellar wind is ellipsoidal, with higher densities on the equator. Followed by the first phase with an integral time of Tw1 = 2.35 × 105 yr, the second stage of the stellar wind lasts for a time of Tw2 = 2.05 × 105 yr with a much larger velocity.

For the Cygnus Loop, the total mass of the ejecta is roughly 4.1 M based on the analysis of the abundance of heavy elements from the observations with ASCA (Miyata & Tsunemi 1999), and our model has Mej = 5.0 M. With an ambient number density of nH = 0.1 cm−3, a subenergetic explosion with a kinetic energy Eej ∼ 0.2 × 1051 erg can lead to a radius consistent with the observations, so we adopt these values in the simulations. As listed in Table 1, the four models with different combinations of vp1, 2, α and ξ are calculated to seek appropriate parameters for the remnant evolved within the cavity.

Table 1.

Parameters for the numerical models with the others ξ = 10, | $\dot{M}_{\rm w} = 10^{-5}$ | M yr−1, T = 104 K, Tw = 4.4 × 105 yr, Tw1 ∼ 2.35 × 105 yr, Tw2 ∼ 2.05 × 105 yr, Mej = 5.0 M, vs = 50 km s−1, nH = 0.1 cm−3 and Eej = 1051 erg.

ParameterModel AModel BModel CModel D
α (°)50505030
ξ1010104
vp1 (km s−1)15151515
vp2 (km s−1)15300150300
ParameterModel AModel BModel CModel D
α (°)50505030
ξ1010104
vp1 (km s−1)15151515
vp2 (km s−1)15300150300
Table 1.

Parameters for the numerical models with the others ξ = 10, | $\dot{M}_{\rm w} = 10^{-5}$ | M yr−1, T = 104 K, Tw = 4.4 × 105 yr, Tw1 ∼ 2.35 × 105 yr, Tw2 ∼ 2.05 × 105 yr, Mej = 5.0 M, vs = 50 km s−1, nH = 0.1 cm−3 and Eej = 1051 erg.

ParameterModel AModel BModel CModel D
α (°)50505030
ξ1010104
vp1 (km s−1)15151515
vp2 (km s−1)15300150300
ParameterModel AModel BModel CModel D
α (°)50505030
ξ1010104
vp1 (km s−1)15151515
vp2 (km s−1)15300150300

The evolution of the stellar wind is indicated with density cuts in Fig. 2 for Model B. During the first stage, a bow shock with an ovoid shape is produced as a result of the interaction between the stellar wind, with α = 50° and vp1 = 15 km s−1, and the ambient medium, which has a velocity of vs = 50 km s−1 towards the −z-direction. The ISM is shocked and is compressed by the shock. The stellar wind flows freely from the star until it encounters the termination shock. Also, a contact discontinuity, which separates the shocked ISM from the shocked stellar wind, is located between the outermost bow shock and the termination shock (Toledo-Roy et al. 2014a). With a wind velocity much smaller than vs, the overall structure of the stellar wind is stable, with laminar flow inside the bow shock (Dgani, van Buren & Noriega-Crespo 1996; Meyer et al. 2015). The stagnation distance derived with equation (1) in Toledo-Roy et al. (2014a) for an isotropic wind is ∼1.2 pc, which is consistent with our calculation. The shock is elongated towards the −z-direction during the evolution, and a low-density cavity is formed between the left and right limbs of the shock. After a time of ∼2 × 105 yr, the border of the shock has arrived at the −z boundary.

Evolution of the stellar wind for Model B. The yz slices of the mass density at x = 0 for different integration times are shown.
Figure 2.

Evolution of the stellar wind for Model B. The yz slices of the mass density at x = 0 for different integration times are shown.

After time Tw1, the wind enters the second phase with a velocity as high as vp2 = 300 km s−1. With such a high wind velocity, much larger than the flow velocity of the ambient medium, the surface of the contact discontinuity is significantly distorted due to Rayleigh–Taylor and/or Kelvin–Helmholtz instabilities (Meyer et al. 2015). The stagnation distance is about | $\sqrt{v_{\rm p2}/v_{\rm p1}} \sim 5$ | times larger than the one of the first phase. As a result, the stellar wind with a higher velocity overpasses the bow shock formed in the first stage, and a new shock with a larger extension is produced in the +z-direction. Finally, the outer shock has a complicated morphology as a result of the interaction of the two-phase stellar wind with the ambient medium.

Assuming that the second phase of the stellar wind lasts for Tw2, the evolution of the remnant is initiated by setting the initial conditions of density, velocity and temperature of the ejecta in the simulations. As indicated in Fig. 3, the ejecta composed of the core and the envelope first expands outwardly, and the density of the core decreases during the expansion. The ejecta expands with a velocity much larger than the sonic speed in the ambient medium, and a shell bordered by a forward and a reverse shock is driven. The shocked medium and the shocked ejecta are separated by a contact discontinuity, where the flow is susceptible to Rayleigh–Taylor instabilities.

Evolution of the SNR for Model B. The yz slices of the mass density at x = 0 for different integration times after the supernova explosion are shown.
Figure 3.

Evolution of the SNR for Model B. The yz slices of the mass density at x = 0 for different integration times after the supernova explosion are shown.

A part of the outward-moving forward shock first interacts with the reverse shock of the bow shock along the direction of the stagnation point, and then it travels into the shocked high-density wind and the shocked ISM. As a result, a bump is generated on the bow shock in the north, and it gradually expands laterally with depressions on both sides. In the south, a bump with an extension of ∼60°, consistent with that observed, is prominent. The extension depends very much on the velocity of the stellar wind during the first phase, and a higher vp1 results in a larger extension of the blowout region with the same spatial velocity of the star. The X-ray morphology of the Cygnus Loop indicates a bright shell on the north-eastern side of the remnant, and it can be explained as the penetration of the forward shock into the bow shock. A bump is seen in the west part of the shell, and the limb of the forward shock from east to south-east is planar with a depression in the middle. These features are consistent with the detected X-ray morphology for the Cygnus Loop remnant.

In Model A, we assumed that the stellar wind has the same velocity distribution vp = 15 km s−1 at all times. As in Model B, a part of the forward shock travels in the undisturbed medium in the +z-direction, whereas the opposite shock evolves in the unshocked wind. Fig. 4 shows that the extension of the bow shock produced by the stellar wind in the +z-direction is much smaller than that in Model B. The resulting profile of the SNR shell in the north-east is more circular than that of Model B. However, the morphology indicated in the slice of the density distribution is inconsistent with that observed, as the observed feature of the limb of the forward shock from east to south is not reproduced in this model. Even in Model C with vp2 = 150 km s−1, the resulting profile in the south-east does not resemble the observations. Therefore, we propose that a stellar wind that has experienced two phases with different velocities is needed to reproduce the detected shock profile for the remnant.

Spatial distributions of the mass density of the stellar wind and the SNR with the integration time of Tw = 4.4 × 105 yr for Models A, C and D, respectively.
Figure 4.

Spatial distributions of the mass density of the stellar wind and the SNR with the integration time of Tw = 4.4 × 105 yr for Models A, C and D, respectively.

An asymmetry of the limbs on the east (left) and west (right) sides of the blowout is prominent in the X-ray morphology of the Cygnus Loop. The west limb is relatively planar, whereas the east limb shows a bump with small curvature radii. The location of the east bump is determined by the inclination angle of the polar axis to the direction of the motion of the star. The panels for Model D in Fig. 4 show the density slices for a stellar wind and the remnant with α = 30° and ξ = 4. With a smaller α, the major axis of the stellar wind is more inclined to the −z-direction, and the east bump is further away from the blowout. With a smaller equator-to-pole density ratio ξ, the matter distribution of the stellar wind is more circular. The bump to the left of the blowout is less extended along the radial direction.

The distributions of the projection of the temperature of the remnant weighted by density along the x-direction (i.e. ∫Tρ dx/∫ρ dx) for the four models are indicated in Fig. 5. The X-rays from the interior of the Cygnus Loop are harder, based on the distribution of the ratio of the hard/soft band, and therefore the inner region is hotter than the edge. For each of the models, the temperature is higher in the inner part of the remnant, which is consistent with the detected hard/soft X-ray band ratio. A loop of higher temperature coincident with the contact discontinuity is more prominent for Model A because the shock profile is less influenced by the ambient matter.

Projected temperatures weighted by density along the x-direction for Models A, B, C and D, respectively.
Figure 5.

Projected temperatures weighted by density along the x-direction for Models A, B, C and D, respectively.

4 SUMMARY AND CONCLUSION

The main goal of this paper is to investigate the formation mechanism of the shock profile of the Cygnus Loop using 3D hydrodynamical simulations. We find that the dynamical properties of the stellar wind originating from the progenitor can be deduced by comparing the mass density distribution with the detected morphology for the remnant. Unlike Toledo-Roy et al. (2014a), in which synthetic X-ray maps are generated to directly compare with the X-ray morphology for Kepler's SNR, we have used the profile of the outermost shock indicated in the X-ray images for the Cygnus Loop to constrain the parameters in the model.

The morphologies seen in the X-ray observations show a blowout in the south and asymmetric limbs on either side of the blowout, a bump in the west and a planar region of the shock from east to south-east. We argue that the progenitor exploded in the circumstellar medium swept by an anisotropic stellar wind, which had experienced two phases with different velocities. To reproduce the asymmetric morphology of the remnant, the progenitor is assumed to move with a velocity of 50 km s−1 relative to the ambient medium, and the stellar wind of the progenitor is anisotropic with an equator-to-pole density ratio of 10. Based on the simulation results, the pole direction of the wind is inclined with respect to the velocity of the star with an angle of ∼50°. Also, the wind has polar velocities of ∼15 and ∼300 km s−1 in the first and second phases, respectively (Model B).

The bow shock of the progenitor continually accumulates the circumstellar medium along the direction of motion of the star to produce a dense bulge, and a low-density cavity is generated in the opposite direction. Afterwards, the dense bulge is swept by the forward shock of the SNR, which can be used to explained the north-eastern bright shell of the Cygnus Loop. Also, the shock in the east has a planar profile and a depression in the middle. In the south, a part of the forward shock expands in the low-density tunnel with wind material, and corresponds to the detected blowout in the south of the remnant. The angular extension of the blowout to the centre of the remnant, the anisotropy of the shell and the relative location and extension of the eastern bump closed by the blowout can be used to constrain the velocity and the anisotropic properties of the stellar wind. After investigating the general evolutions of the wind and the remnant, we have obtained a morphology of the shock comparable to that seen in the X-ray band using appropriate parameters for the SNR and the stellar wind. We have not explored the whole parameter space in this paper because the 3D hydrodynamical simulations are very time-consuming.

We thank an anonymous referee for insightful comments that have improved this work. JF acknowledges the support of the National Natural Science Foundation of China (NSFC; 11563009), the Open Research Programme of the Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences (OP201402), the Yunnan Applied Basic Research Projects (2016FB001) and the Top Talents Programme of Yunnan Province, China (2015HA030). LZ acknowledges support from NSFC 11433004. HY is partially supported by the Yunnan Applied Basic Research Projects (2016FD105) and IRTSTYN.

REFERENCES

Aschenbach
B.
Leahy
D. A.
,
1999
,
A&A
,
341
,
602

Blair
W. P.
Sankrit
R.
Raymond
J. C.
,
2005
,
AJ
,
129
,
2268

Bocchino
F.
Orlando
S.
Miceli
M.
Petruk
O.
,
2011
,
A&A
,
531
,
A129

Dalgarno
A.
McCray
R. A.
,
1972
,
ARA&A
,
10
,
375

Dgani
R.
van Buren
D.
Noriega-Crespo
A.
,
1996
,
ApJ
,
461
,
927

Fang
J.
Zhang
L.
,
2012
,
MNRAS
,
424
,
2811

Fang
J.
Yu
H.
Zhang
L.
,
2014
,
MNRAS
,
445
,
2484

Fraschetti
F.
Teyssier
R.
Ballet
J.
Decourchelle
A.
,
2010
,
A&A
,
515
,
A104

Hester
J. J.
Cox
D. P.
,
1986
,
ApJ
,
300
,
675

Jun
B.-I.
Norman
M. L.
,
1996
,
ApJ
,
465
,
800

Levenson
N. A.
et al. ,
1997
,
ApJ
,
484
,
304

Levenson
N. A.
Graham
J. R.
Keller
L. D.
Richter
M. J.
,
1998
,
ApJS
,
118
,
541

Levenson
N. A.
Graham
J. R.
Snowden
S. L.
,
1999
,
ApJ
,
526
,
874

McEntaffer
R. L.
Brantseg
T.
,
2011
,
ApJ
,
730
,
99

Meyer
D. M.-A.
Langer
N.
Mackey
J.
Velázquez
P. F.
Gusdorf
A.
,
2015
,
MNRAS
,
450
,
3080

Mignone
A.
Bodo
G.
Massaglia
S.
Matsakos
T.
Tesileanu
O.
Zanni
C.
Ferrari
A.
,
2007
,
ApJS
,
170
,
228

Mignone
A.
Zanni
C.
Tzeferacos
P.
van Straalen
B.
Colella
P.
Bodo
G.
,
2012
,
ApJS
,
198
,
7

Miyata
E.
Tsunemi
H.
,
1999
,
ApJ
,
525
,
305

Orlando
S.
Bocchino
F.
Reale
F.
Peres
G.
Petruk
O.
,
2007
,
A&A
,
470
,
927

Orlando
S.
Petruk
O.
Bocchino
F.
Miceli
M.
,
2011
,
A&A
,
526
,
A129

Orlando
S.
Miceli
M.
Pumo
M. L.
Bocchino
F.
,
2016
,
ApJ
,
822
,
22

Raga
A. C.
Cantó
J.
De Colle
F.
Esquivel
A.
Kajdic
P.
Rodríguez-González
A.
Velázquez
P. F.
,
2008
,
ApJ
,
680
,
L45

Schneiter
E. M.
Velázquez
P. F.
Reynoso
E. M.
Esquivel
A.
De Colle
F.
,
2015
,
MNRAS
,
449
,
88

Toledo-Roy
J. C.
Velázquez
P. F.
de Colle
F.
González
R. F.
Reynoso
E. M.
Kurtz
S. E.
Reyes-Iturbide
J.
,
2009
,
MNRAS
,
395
,
351

Toledo-Roy
J. C.
Esquivel
A.
Velázquez
P. F.
Reynoso
E. M.
,
2014a
,
MNRAS
,
442
,
229

Toledo-Roy
J. C.
Velázquez
P. F.
Esquivel
A.
Giacani
E.
,
2014b
,
MNRAS
,
437
,
898

Tsebrenko
D.
Soker
N.
,
2015
,
MNRAS
,
450
,
1399

Tsunemi
H.
Katsuda
S.
Norbert
N.
Miller
E. D.
,
2007
,
ApJ
,
671
,
1717

Uchida
H.
Tsunemi
H.
Katsuda
S.
Kimura
M.
,
2008
,
ApJ
,
688
,
1103

Yu
H.
Fang
J.
Zhang
P. F.
Zhang
Li.
,
2015
,
A&A
,
579
,
A35

Zhou
X.
Bocchino
F.
Miceli
M.
Orlando
S.
Chen
Y.
,
2010
,
MNRAS
,
406
,
223