Modelling colliding wind binaries with RAMSES, extension to special relativity
Astrid Lamberts, Sebastien Fromang, Guillaume Dubus, Geoffroy Lesur
***Volume Title**ASP Conference Series, Vol. **Volume Number****Author** c (cid:13) **Copyright Year** Astronomical Society of the Pacific Modelling colliding wind binaries with RAMSES, extension tospecial relativity
Astrid Lamberts, S´ebastien Fromang, Guillaume Dubus, and Geo ff royLesur UJF-Grenoble 1 / CNRS-INSU, Institut de Plan´etologie et d’Astrophysique deGrenoble (IPAG)UMR 5274, Grenoble, F-38041 Laboratoire AIMCEA / DSM-CNRS-Universit´e Paris Diderot, IRFU / Service d’AstrophysiqueCEA-Saclay F-91191 Gif-sur-Yvette, France
Abstract.
We present simulations of colliding supersonic stellar winds with RAM-SES. The collision results in a double shock structure that is subject to di ff erent instabil-ities.Properly modelling these instabilities requires a high enough resolution. At largescale, orbital motion is expected to turn the shocked zone into a spiral but we find thatin some configurations the Kelvin-Helmholtz instability (KHI) may disrupt the spiral.A colliding wind structure is also expected in gamma-ray binaries composed of a mas-sive star and a young pulsar that emits a highly relativistic wind. Numerical simulationsare necessary to understand the geometry of such systems and should take into accountthe relativistic nature of the pulsar wind. We implemented a second order Godunovmethod to solve the equations of relativistic hydrodynamics (RHD) with RAMSES, in-cluding the possibility for Adaptive Mesh Refinement (AMR). After a brief overview ofour numerical implementation, we will present preliminary simulations of gamma-raybinaries.
1. Introduction
Massive stars posses highly supersonic winds. In binary systems, the interaction oftwo winds creates a double shock structure which geometry depends on the momentumflux ratio of the winds. The winds are separated by a contact discontinuity. Numericalsimulations show that several instabilities may arise in the colliding wind region (Pit-tard 2009). For isothermal winds, the non-linear thin-shell instability (Vishniac 1994)dominates the collision region while the KHI can strongly perturb the shocked regionbetween adiabatic winds (Lamberts et al. 2011).Up to now, high-resolution simulations have focused on the region close to thebinary. At larger scale orbital motion turns the shocked structure into a spiral. Matteris expected to be ballistic, with a distinction between both shocked spiral arms (Parkinet al. 2011). The exact structure of the spiral is still to be studied. A spiral structureis observed in WR 104, a binary composed of a Wolf-Rayet (WR) and an early-typestar. Owing to dust emission, it displays a spiral in infrared up to a few hundred timesthe binary separation (Tuthill et al. 2008). The formation of dust in such systems isstill poorly understood, and is likely to be facilitated in the colliding wind region. We1 a r X i v : . [ a s t r o - ph . H E ] D ec Lamberts, Fromang, Dubus and Lesur perform simulations of this system to determine the structure of the colliding windregion and put constrains on dust formation. γ -ray binaries are binaries composed of a massive star and a young pulsar (seee.g. Abdo & Fermi Collaboration (2009)) with a tenuous, highly relativistic wind. Thecollision between the pulsar wind and the wind from the companion star, which re-sults in γ -ray emission, is expected to present similarities with the collision betweenstellar winds (Dubus 2006). Still, the relativistic nature of the pulsar wind is likely tomodify the dynamics and emission properties of large scale colliding wind region thatcan be observed in radio (see e.g. Mold´on et al. (2011)). To investigate the impactof the relativistic nature of the pulsar wind, we developped a relativistic extension tothe hydrodynamical code RAMSES (Teyssier 2002). We explain our numerical imple-mentation, focusing on the AMR and present some preliminary simulations of γ -raybinaries.
2. Large scale structure of colliding wind binaries
We use the hydrodynamical code RAMSES for our simulations. It is a second orderGodunov method, that solves the equations of hydrodynamics. We focus on the adia-batic limit. We use AMR to increase resolution at the shocks and properly model theinstabilities while simulating the large scale structure at reasonable computational cost.The winds are generated following the method by Lemaster et al. (2007). We include apassive scalar to distinguish the winds and determine mixing.
Figure 1. Density (upper row) and mixing (bottom row) for 2D simulations withincreasing velocity ratios and decreasing density ratios (from left to right). On theleft panel both winds are identical. The length scale is the binary separation. Themaximum resolution is set by l max = amberts, Fromang, Dubus and Lesur M v = ˙ M v ) but increasing velocity ratios (and decreasing mass loss rate ratios). When bothwinds are identical (left panel), the spiral arms have the same step and same width.When the winds have dierent speeds, the spiral arm propagating into the faster, lowerdensity wind expands, while the arm propagating into the denser and slower wind getscompressed. Mixing is more important in the wider arm (Lamberts et al. 2012). Whenincreasing the velocity dierence even more ( v = v = v × P where P is the orbital period ofthe system and v the speed of the dominating wind, with significant corrections fromthe weaker wind in some cases.We perform a detailed study of WR 104. Our goal is to determine whether a hy-drodynamic model with adiabatic winds can reproduce the observed spiral structure inWR 104. The WR wind is very hostile to dust formation because of its low density, hightemperature and absence of hydrogen and the colliding wind region is likely to facilitateit (Walder & Folini 2003). In WR 104 the wind from the WR strongly dominates thewind from the early-type star and the shocks form very close to the latter. A very highresolution is thus necessary close to the binary and is incompatible with a large scale3D simulation. We thus perform a small scale 3D simulation and a complementarylarge 2D simulation. The 2D simulation confirms the Archimedean spiral and indicatesthe highest densities are reached at the edges of the spiral, suggesting they could bethe location for dust formation (Lamberts et al. 2012). The results from the 3D sim-ulation are given on Fig.2. The temperature map suggests temperature could becomelow enough to enable dust condensation ( ≤ (cid:39) cm − ). An adiabatic model is thus unable to account for dust formation in WR104. Cooling is probably present, and increases the compression ratio of the shocks.It might also trigger the non-linear thin shell instability, which results in small highdensity zones. Figure 2. Density (g cm − ), mixing and temperature (K) in the orbital plane ofthe 3D simulation of WR 104. The length scale is the binary separation. Lamberts, Fromang, Dubus and Lesur
3. Extension of RAMSES to relativistic hydrodynamics, application to γ -ray bi-naries The equations of RHD can be written as a system of conservation equations ( c ≡
1) : ∂ D ∂ t + ∂ ( Dv j ) ∂ x j = ∂ M i ∂ t + ∂ ( M i v j + P δ ij ) ∂ x j = ∂ E ∂ t + ∂ ( E + P ) v j ∂ x j = DM i E = Γ ρ Γ ρ hv i Γ ρ h − P (1)where D is the mass density, M the momentum density and E the energy density inthe frame of the laboratory. The subscripts i , j stand for the dimensions, δ i , j is theKronecker symbol. h is the specific enthalpy, ρ is the proper mass density, v i is the fluidthree-velocity, P is the gas pressure and γ the adiabatic index. The Lorentz factor Γ isgiven by Γ = √ − v (2)These equations have a similar structure to the equations of hydrodynamics but aremore complex to solve because strongly coupled to each other by the Lorentz factorand the enthalpy. An additional numerical constraint arises from the the fact that thevelocity must remain subluminal. The similarity with the equations of hydrodynam-ics allows us to closely follow the algorithm implemented in RAMSES, performinglocalised changes.A first di ffi culty arises when passing from the conservative variables ( D , M , E ) T to the primitive variables ( ρ, v , P ) T . There is no explicit solution, and a root find-ing algorithm is necessary to recover the primitive variables. Care has to be taken toavoid numerical problems in the ultrarelativistic and non-relativistic limits (Mignone& McKinney 2007). Second-order precision is implemented in RAMSES following aMUSCL-Hancock method. The reconstruction of the primitive variables in space is nota ff ected by special relativity but the prediction in time requires the determination ofthe Jacobian matrix of the system given by Eq. 1. In case of a superluminal velocityin the reconstructed state, the code switches to a first order scheme. The relativisticsummation of velocities changes the determination of the wavespeeds and the timestep.In RAMSES, AMR is implemented following a tree-based method where par-ent cells are refined in child cells in a recursive structure (Kravtsov et al. 1997). Thechild cells are gathered together in octs of size 2ndim. Interpolation from l to l + D > , E > D + M ). Hydrodynamical updates are only performed at thehighest level of refinement, and variables are computed at lower levels by averagingvalues over the child oct. Although each cell from the child oct at level l satisfies( D > , E > D + M ), it is not necessarily true on average, over the whole oct. Thismay lead to non-physical states at level l −
1. Therefore, we choose to perform therestriction on the internal energy rather than the total energy. This method is currentlyimplemented in RAMSES.Fig. 3 shows two examples of numerical tests we performed with our new rel-ativistic code. The left panel shows the density, transverse and parallel velocity and amberts, Fromang, Dubus and Lesur Figure 3. Left panel : Sod test ( t = . Γ max = . t = , , pressure in a one-dimensional Sod test. Contrary to the classical case, transverse ve-locities do impact the structure of shocks in RHD. The setup is taken from (Ryu et al.2006) and constitutes a very stringent test due to the high Lorentz factor ( Γ max = γ -ray binaries. We focus on the small scale structure ofthe interaction between a stellar wind and a pulsar wind. The goal is to understandthe impact of relativistic e ff ects both on the structure and stability of the interactionregion.We neglect orbital motion and focus on winds with equal moment fluxes. Weperform preliminary simulations with various values of the momentum flux ratio andpulsar wind velocity. They prepare a large scale simulation to determine whether astable structure is possible (Bosch-Ramon et al. 2012). Fig. 4 shows the density mapfor a simulation with the speed of the pulsar wind v p = .
5. The KHI develops in asimilar fashion than in the classical case. The right panel shows the positions of thediscontinuities for simulations with increasing values for the pulsar wind speed. Thehigher the value, the more the shocks are bent towards the star. This is a relativistice ff ect due to the impact of transverse velocities on the structure of shocks.
4. Conclusions and perspectives
We performed high resolution simulations of colliding wind binaries at a spatial scalenever reached before. We showed that the KHI may destroy the expected large scalestructure. Simulations of WR 104 match well with the observed structure and indicate
Lamberts, Fromang, Dubus and Lesur
Figure 4. Left panel : density map of a simulation with equal momentum flux,with v p = ,
5, the star is on the left, the pulsar on the right. Right panel : position ofboth shocks and the contact discontinuity, in simulations with di ff erent values for thevelocity of the pulsar wind. We have v p = .
01 (thin dotted line), 0.1 (thick dottedlined), 0.5 (thick dashed line) and 0.9 (solid line). cooling has to be taken into account to allow dust formation in this system. To model γ -ray binaries, we extended RAMSES to relativistic hydrodynamics. Preliminary sim-ulations of γ -ray binaries confirm a similar structure to stellar binaries. The relativisticextension of RAMSES allows the use of AMR and is suited for the study of gamma-raybursts, relativistic jets or pulsar wind nebulae. It will be part of the next public release. Acknowledgments.
AL and GD are supported by the European Community viacontract ERC-StG-200911. Calculations have been performed at CEA on the DAPHPCcluster and using HPC resources from GENCI- [CINES] (Grant 2011046391)
References
Abdo, A., & Fermi Collaboration 2009, ApJ, 706, L56.
Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59.
Del Zanna, L., & Bucciantini, N. 2002, A&A, 390, 1177. arXiv:astro-ph/0205290
Dubus, G. 2006, A&A, 456, 801. arXiv:astro-ph/0605287
Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73. arXiv:astro-ph/9701195
Lamberts, A., Dubus, G., Lesur, G., & Fromang, S. 2012, A&ALamberts, A., Fromang, S., & Dubus, G. 2011, MNRAS, 418, 2618Lemaster, M. N., Stone, J. M., & Gardiner, T. A. 2007, ApJ, 662, 582. arXiv:astro-ph/0702425
Mignone, A., & McKinney, J. C. 2007, MNRAS, 378, 1118.
Mold´on, J., Johnston, S., Rib´o, M., Paredes, J. M., & Deller, A. T. 2011, ApJ, 732, L10 + . Parkin, E. R., Pittard, J. M., Corcoran, M. F., & Hamaguchi, K. 2011, ApJ, 726, 105.
Pittard, J. M. 2009, MNRAS, 396, 1743.
Ryu, D., Chattopadhyay, I., & Choi, E. 2006, ApJS, 166, 410. arXiv:astro-ph/0605550
Teyssier, R. 2002, A&A, 385, 337. arXiv:astro-ph/0111367
Tuthill, P. G., Monnier, J. D., Lawrance, N., Danchi, W. C., Owocki, S. P., & Gayley, K. G.2008, ApJ, 675, 698.
Vishniac, E. T. 1994, ApJ, 428, 186. arXiv:astro-ph/9306025arXiv:astro-ph/9306025